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ABSTRACT 


Exact solutions for the flow of a compressible laminar boundary 
layer over a flat plate with homogeneous suction or injection are 
considered in this report 
given in detail, and it is shown that the solutions by the momen 
tum method given in references 1, 2, and 9 compare favorably 
with these results 

Asymptotic solutions are also developed which are valid for 


the flow far downstream from the leading edge of the plate or for 


flows over arbitrary bodies with large suction It is shown that 
the exact solutions for the flow over a flat plate with homogeneous 
suction approach the asymptotic solution downstream of the 
leading edge. Moreover, for practical purposes, this occurs at 
arather short distance from the leading edge. The asymptotic 
solutions are given in closed form for several cases; they may also 
be used as a first approximation for the above flows. It is also 
shown that the method of Lighthill'? for the approximate cal 
culation of heat transfer at the wall 1s exact for the asymptotic 


solution 


SYMBOLS AND NOTATION 


Pi (Cur NV 
constant defined by u* 
specific heat at constant pressure 
Cpl + (u?/2 
Col 
coefficient of thermal conductivity 
Mach Number 
(y — 1)/2]A7,? 
= Sutherland constant 
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Characteristics of these flows are 


absolute temperature 

ie § 

free-stream velocity in x direction 

velocity components in x and y directions, respectively 

coordinates parallel and perpendicular to surface, re 
spectively 

Nusselt Number 

Prandtl Number, (uc,,/k 

distance parameter, | — 

Vv 2Cy,x/ l 


distance parameter, (—7/ Ul")? (Ux/1 


distance parameter, | —(p Cu 


ratio of specific heat at constant pressure to specific 
heat at constant volume 
kinematic viscosity (u/p 
mass density of fluid 
coeflicient of viscosit\ 
shear stress, also dependent variable in Eq 
stream function 
p” p/ pi 
“ M/E 
Subscripts 
& conditions outside of the boundary layer 
0) = conditions at the wall or surface of the body 
Superscript 


nm = exponent occurring in p* a ie 


INTRODUCTION 


A STUDY OF THE FLOW of a compressible boundary 
layer over a flat plate with a small uniform normal 
velocity at the plate is presented in this report. This 
solution might be found useful in several applications. 
On the one hand, the use of suction to control the 
boundary layer at high speeds requires the accurate 
knowledge of the velocity and temperature profiles 
On the other hand, the possibility of the use of injection 
at the wall to cool a body in a high-temperature gas 
demands a knowledge of the heat-transfer propertics 


in such a combination. Several reports have dealt 
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Ux,y) 


GTM TTT 


Vo\x) = CONSTANT 


Fic. 1 Explanatory sketch for flat plate in longitudinal flow 
with constant suction 





with solutions of the above problems. For the case of 
suction, which this report will deal with for the most 
part, Lew! has solved the flow of a compressible bound- 
ary layer over a flat plate with homogeneous suction 
by a modified Karman momentum method. In the 
case of homogeneous injection, the problem has been 
solved by Lew? utilizing the modified Karman momen- 
tum method while Yuan* solved the same problem by 
applying the Karman-Polhausen method directly as 
in the incompressible flow. In this report, the prob- 
lem is explored and exact solutions are calculated for 
the cases of homogeneous suction and_ injection. 
Asymptotic solutions for the suction case are derived 
in closed form for the case of arbitrary heat transfer at 
the boundary. These asymptotic solutions represent 
the flow at an infinite distance from the leading edge 
in the case of homogeneous suction and the flows for 
extremely large suction which are compatible with the 
boundary-layer assumptions. Moreover, they showed 
the effects of the various parameters easily due to their 
simplicity of solution. The asymptotic solutions also 
indicate the validity of an approximate method for the 
calculation of heat transfer at the boundary. 

This work was done with the thought that an exact 
solution might be found useful for comparison with the 
approximate solutions and especially for stability 
calculations. Extensions to the cases of an arbi- 
trary pressure gradient and arbitrary suction or injec 
tion distribution could be made. 


EQUATIONS OF MOTION 


An x-y-Cartesian coordinate system with the origin 
at the leading edge of a semi-infinite plate is utilized 
here (see Fig. 1). The equations of the compressible 
boundary layer are the momentum, mass, and energy 
conservation laws, and the equation of state of a per 
fect gas. They are, respectively: 


pt [u(Ou/Ox)] + [v(Ou dy) ]{ = 
(O/Oy) [u(Ou/Oy)] (1) 


[O(pu)/Ox] + [O(pv) Ov] = O (2) 

( oT in = o (: =) 4 (*) (3 

Ip iu v = 2 » 
ee ox | Oy) dy\" Oy Noy 

Pp -_ pRT (4) 


A linear viscosity temperature relation is utilized here 


of the form 


MMi Clit / Fs 


as given in reference 4. The constant C 1s determined 
by the more accurate Sutherland relation for the ten 
perature dependence of viscosity. The use of this 
relation [Eq. (5)] has been discussed in reference 

{n addition, the Prandtl Number is taken as unity in 


the numerical calculations so that some former results 
can be used. However, this assumption is not necessary 
in the analysis. The extension to arbitrary Prandtl 
Number is given in Appendix II of reference 13. For 
a Prandtl Number of unity, the familiar Crocco i1 
tegral of the energy equation could be used as shown j 
reference 6. The momentum and energy equations 
are satisfied by a polynomial relation between the 


temperature and the velocity of the form 


In addition, if there is zero heat transfer at the wall 


then (07° Oy), = O, and since the free-stream total 
temperature ratio is given by 
10/71 = 1 + ([(9 2 LA; 


Eq. (6) reduces to 
1 /fy = (16/70) (y — 1)/2).Mir(a/ Ll 7 


The boundary conditions for this problem are 


u(x, 0) = O 
nix, o) = { 
u(0O, y) = l 
\ 
V(x, O) = vo(X 
T(x, 0) = Ty = constant 
I(x, o) = T, = constant 


The above equations can be more convenicitly used 
in the (x, yw) 
equation |Eq. (1) ] becomes 


coordinates. Thus, the momentum 


Ou/Ox = (1/ pi") (0 OW) [upu(du Oy 9 

where the stream function y is defined by 
u = (p/p) (Ow, Oy / ” 

= (p:p) (Ow /Ox) \ 


The boundary conditions on the velocity components 


become 


ulx, f(x) | ri 
u(x, 2 l 1 
u(O, l | 
and 
* 
f(x) - | Vo(x)dx it 4 0 
7) Pi 
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AMINAR COMPRESSIBLE BOUNDARY LAYER OVER A FLAT PLATE 


x 


The boundaries of integration are curvilinear as 
can be seen from Eqs. 1] and the following trans 
formats leads to straight line boundaries that are 

ore ible to analysis. Let 

2) y Aw 
> \ |? 
/ \ { Pe 
n Rae 
J \ t [ dx | 
where is as defined by eos. (11) and the prime de 


notes differentiation. Thus, Eq. (9) becomes 


Olt f(x)? Ou 


with the boundary conditions 


u(n, O () 
u (n, l 14 
u(O, d [ | 


For the case of constant v(x), Eq. (13) and boundary 


Eqs. (14)] can be further transformed to 


conditions 


the more convenient forms of 


It) 


The new variables are defined by 


r V —A.o/V (A,2/Bi)n | 


IS 
V (A,? By)n 
whert 
| Pui Cu 
By = pi/(Cu) V1 + (p0?/ pr?) 0?/U 
Eq. (15) and the boundary conditions |Eqs. (16)] are 


valid for the flow over a flat plate with uniform suction 
It is noted that if suction is considered 
Corre 


takes on negative values if injection is 


OT injection 
the constant A, is negative and a is positive. 
spondingly, 
In these relations the viscosity relation as 
In this form, Eqs. 


considers 
has been used. 
16) are independent of the magnitude of the 


given by Eq rs) 
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suction or injection velocities at the wall and of the 
physical properties of the fluid Moreover, the equa 
tions have the same form as in the incompressible case 
so that all results of the incompressible case can be 
used 

The conversion to the physical plane can be easily 


acce miplished by the relation 


rua (7 [ 1 7) TUT F 
2 - | 2o lor constant o (19 
Cp zy JO / V 


where ¢€ (—w Ll lx vy Phe second of Eqs 
19) is found from Eq. (10) by integration in the y 
direction for constant x. 

At the point ¢ Q for either suction or injection 
there is a singularity of the differential equation Eq 
(15 However, it can be shown that the velocity 
profile at this point 1s a Blasius profile If an expan 
sion of I’(o, 7) 1 a series in the neighborhood of o 0 


is inade of the form 


Vio, 7 | T al yr T) + 20) 
then Eq. (15) leads to 
VS Von” + Vo'[r VVo/tT 0 2] 
with the insertion of the expansion, Eq. (20 The 
boundary conditions for I'y(7) are 
1,(0 0, V } 22 


It is seen that Eqs. (21) and (22) are the Blasius system 
in the o, 7 coordinates, so it follows that the Blasius 
profile exists at o 0 for both suction and injection 
Iglisch has shown this for the suction case’ for incom 
pressible flow. 

It may be noted that the solutions of the compressible 
boundary-layer equations with suction and zero pres 
sure gradient can be deduced from the incompressible 
one when yu is proportional to 7° by using the trans 
formation given by the senior author in references | 
and 9 as a generalization of Dorodnitzyn’s work for 
zero heat transfer. This transformation is a special 
case of Illingworth and Stewartson’'s results given in 
reference 14 and 15 on the equivalence of compressible 
and incompressible boundary laver for the cases with 


non-zero pressure gradient 


EVALUATION OF THE CHARACTERISTICS OF THE FLOW 
The characteristics of the boundary layer are con 
veniently given in terms of the variable y in the form of 
y wy) /(Cv 23 
The first and second derivatives of the velocity com 
ponent w are then given by (a 
1 Ou T, 7; OV 
l Ov" ST, Tor Or 


and (b 
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Fic. 2. Skin friction vs. distance parameter £ for suction and 
injection 


1 Ou ~v SERIES bees () 

U dy*? ~ 1602r\To/ \T /L 7dr Or \T 
| Cr) ol =~) | 
PAT /\" de Or 


Thus the skin friction at the wall given by the ex- 


pression 


T) = Mo(OU/ OY) (25) 

iS 
T) = — pot) U(do/4e) (26) 
where d. = (1/2) (07?V/0Or?),=9 and o is related to é 
by the first of Eqs. (19). Eq. (26) is shown in Fig. 2. 


(Cf. Table 1.) The displacement thickness defined by 


6* = [ (: «= * | dy (27) 
J0 pi 


becomes a function of o in the form of 


~ = ae o [ Ie = Vv V(o, a] x 
( Vi T; 70 T; r 2 


rdr 


Vv V(o, T) 


(28) 


The integration with respect to 7 is performed with o 
held constant. Lastly, the momentum thickness de- 


6 2 “(1 ) dy (29) 
v0 pi l l : 


is given in the form of 


fined by 


V0 T 


. ae [ ( - : V(o, 2} rdr (30) 
Cr, i; 0 D) 


CasE OF HOMOGENEOUS SUCTION 
The solution of Eq. (15) with the boundary condi- 
tions [Eqs. (16)] for the suction case (positive a) has 
been calculated in detail by Iglisch® for incompressible 
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flow by a method of profile continuation using Adam's 
method of numerical integration for a fixed value oj 
o. The integration involves replacing the differential 
term in the o direction by a finite difference approxi- 
mation and linearizing the equation by replacing all 
terms of ~/ V by values at the previous value of ¢ 
The resulting ordinary differential equation is then 
integrated by Adam’s method. From Iglisch’s results 
the values of V’(o, 7) for the case of suction have been 
calculated. The suction case, of course, corresponds 
to ao > 0. These values are tabulated in reference 13 
From these values of J’(¢, 7), the various character 
istics in the compressible case can be calculated by 
one integration for y* in the form of the first-order 
differential equation 


dy*/dr = Q(r) (31 


with the initial condition of y* = 0 at 7 = 0. Here 
Q(7r) is given by the second of Eqs. (19) in the form 


O(r) = (270/71?) o[(T7)/V Ve, r)] 


where o is held constant. The temperature 7° is re 
lated to Il’ by Eqs. (6) or (7) and the first of Eqs. (17 

The velocity profiles at .1/,; = 4 for several values of 
£ have been calculated for zero heat transfer at the wall 
and are given in Fig. 3. For the purpose of comparison 
with reference 9, which uses a modified momentum 
method, one velocity profile was calculated for 1, = 
2, 10/11 = Zaté 0.5 for C 1. This calculation is 
shown in Fig. 4 where the result of the momentum 
method (reference 9, for 1 = 1) is also shown. The 
comparison is extremely good indicating that the 
momentum method is accurate as far as the velocity 
profiles are concerned. A comparison of the first de- 
rivatives of the velocity profiles is given in Fig. 5, and 
again the agreement between the exact method and 
momentum method is excellent. As to be expected, 
the second derivatives of the profiles using the momen- 
tum method do not agree well with the exact method 
given here. Values of the second derivatives for the 
momentum method differ by about 20 per cent at the 


TABLE | 
d. Coefficients for Various o 


o ds» 

0.2 1.412 
0.1 1.633 
0 1.8738 
O.1 2.129 
0.2 2.389 
0.3 2.659 
0.4 2 936 
0.5 3.224 
0.6 3.521 
0.7 3.825 
O.8 4.134 
0.9 +. 453 
1.0 4.776 
1.2 5.448 
1.4 6.126 
:.6 6.835 
2.0 § 288 
2.4 9.782 
3.2 12.915 
4.0 16.000 
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Fic. 3. Suction velocity profile development . ; 
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Fic. 4. Comparison of suction velocity profile with reference 9 
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Fic. 7. Comparison of injection velocity profile with reference 2 


wall. A comparison with reference 9 of the first de- 
rivative at the wall as a function of o for 1/4; = O is 


shown in Fig. 6. 


CASE OF HOMOGENEOUS INJECTION 


For the case of homogeneous injection, Eq. (15) 
with boundary conditions [Eqs. (16)] is integrated for 
The same method as described 
Detailed calcula- 


negative o from a = 0. 
under the suction section is used. 
tions are presented from ¢ = 0 to o = —O.2 in refer- 
ence 13. The comparisons of the characteristics of the 
flow with the momentum method? are shown in Figs. 
7 and 8. Again the velocity profiles and its first de- 
rivatives are in good agreement even though the com- 
parison is made at a point extremely close to the leading 
edge of the plate where characteristics of the flow vary 
rapidly. The second derivatives do not compare 
well; for example, they differ from each other by about 
17 per cent at the wall. The Nusselt Number char- 
acterizing the heat transfer at the wall defined by 


N, = —(OT/Oy)o [x/(T1 — To)] 
appropriate to this analysis is 


N, = = — 
4aC(Ty i) [1 _ (To, T;)] 


9) 


At € = 0.005, T/T; = 1/2, M, = 2, C = 1, Eg. (32) 
gives a value of 


Ny = —0.09177(U/u) 
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as compared to the momentum method? value of NV, = 
—0.10048(U/u). The difference is 8.7 per cent for 
this value of £, which is very near the leading edg 
The skin-friction intensity is also shown in Fig. 2, indj. 
cating a smaller value than the suction case, as is to 
be expected. Moreover, this skin-friction intensity 
will approach zero, indicating a possible separation of 
the flow. Details of the calculations are given jp 
Appendix I of reference 13. 


ASYMPTOTIC SOLUTIONS 


The asymptotic characteristics of the flow in the case 
of suction will be considered in this section. These 
asymptotic properties will hold for large values of ; 
or large values of the suction velocity at the boundary 
compatible with the boundary-layer assumptions 
The asymptotic results are important since the flow 
for suction will rapidly approach this state a short 
length from the leading edge of the plate. Moreover, 
these equations are some of the relatively few in the 
theory of compressible laminar boundary layer which 
can be integrated in closed form. The following case 
is first considered with the Prandtl Number equal to 
unity. If the conditions far downstream from the 
leading edge of the plate are considered, then all proper- 
ties of the flow are independent of x. Hence Eqs 


(1), (2), and (3) become 


pu(du,dy) = (d/dy) [u(du/dy) | 30 
pv = constant (34 
pu(dT'/dy) = (1/P,) (d/dy) [u(dT/dy)| + 
(u/cy) (du/dy)? (35 
for P, = constant. The boundary conditions are 
y = 0: u = 0, T = T, = constant / 
(36 
y> mo: u=U,T=T, \ 


When the Prandtl Number is unity the temperature 
is related to the velocity by Eq. (6). In this analysis, 
the viscosity can be taken as a function of the temper- 
ature in the form 


~ 


u* = CT*" (37 
Thus for m = 1, Eq. (37) coincides with Eq. (5). For 
c = 1 the usual power relation is obtained for which 


is taken usually at 0.7 to 0.8 for air. 
Eqs. (33) to (37) can be integrated to the form 


me [ “(A + Bt + Ne?)" di (38) 
: pv J0 tk = 2) 
where 
A = 1/T; 
B = {1 — (%/T) + [(v — 1)/2] MV} 
N = —[(y7 — 1)/2])M?? 


Here ¢ is a dummy variable of integration. For » = 1, 


Eq. (38) integrates to 
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c) = (To/T1) (CN /2) (u/U)? + 
(B + N) (u/U) + log{1 — (u/U)]f (39 


The boundary conditions [Eqs. (36)] have been satis- 
fed in Eqs. (38) and (39). It is noted that Eq. (38 
an be simply integrated numerically for nonintegral 
values of nm. Equation (39) reduces to the incompres- 
sible asymptotic profile of Griffith” if 4 = C = | 
nd M, = N =0. 

For Prandtl Number not equal to unity, an extension 
{the solution by Young,'! who considered the problem 
with zero heat transfer at the wall, to the case of non- 
zero heat transfer at the wall is possible. Equations 
23) to (35) are first transformed from the (y) coordi- 


nate to the (uw) coordinate; they become 
u(du/dy) = —pw(U — u) (40) 
1 — P,) (di/du) — (U — u) X 
1 [(d22)/(du2)]} + P,f =O (Al 
with the boundary conditions of 


ty aty =0 | 


II 


u=Q, 1 
(42) 


u=U, t= as y—> o| 


It can be shown that the solution to Eq. (41) and satis- 
fying conditions (42) for arbitrary Prandtl Number is 


P, a 
_ _ ou —— 3" {2 4 
=: 9 (: *) (1 lo 7 9 ) 
PU? ! =) re ( an 
(1 =— — — l- = (453 
2(2 — P,)L\ l P, I 


The velocity w is related to y by 


u/U = 1 — exp(¢) (44) 
and 
” i n 1 : 
c= - : — (45) 
: J0 1 és 
where » = —(uy)/m. The viscosity relation (37) is 
used here. Substitution of Eq. (44) into Eq. (43) leads 
to 
1— T*)/N = (B/N) exp (Pt) + 
[P,/(2 — P,)|]F(¢) (46) 
where 


F(¢) = (2/P,) exp (Pf) — exp (2¢) 


For the case where the viscosity varies linearly with the 
temperature so that m = 1, Eq. (45) can be integrated 


in closed form and is 


” Z we i 
ee ( a 
P 


‘) = 
(47) 


In fact, Eq. (45) can be integrated in closed form for n 
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Fic. 8. Comparison of first derivative for injection profile with 


reference 2 


equal to an integer. Equation (47) reduces to Eq. 
(39) for the case of P, = 1. 


perature gives an indication of the heat transfer at the 


The gradient of the tem 


wall; thus, 

(dT* dy)o = [(vB) (1;C) | (7; ig? Pr, (48) 
using Eqs. (45) and (46). The skin-friction intensity 
is particularly simple and has the same form as in the 


incompressible case of 
T = —pol'm (49) 


The temperature profiles for various conditions with 
n = 1 are shown in Fig. 9. The effect of Prandtl 
Number on the temperature profiles is pronounced, as 
seen from these curves that were calculated for P, = 
1 and P, = 0.72. In fact, the temperature gradient 
[Eq. (48)] at the wall is directly proportional to the 
value of the Prandtl Number. Its value for P, = 1 
is 40 per cent in error if P, = 0.72 is assumed to be 
correct. This is another indication that the heat- 
transfer characteristics should be based on correct 
values of P,. 

The effect of Prandtl Number on the velocity pro- 
files is negligible as can be readily seen from Fig. 10. 
The change of the suction velocity profiles from the 
leading edge to the asymptotic profile occurs rapidly 


(Fig. 3). For example, at a value of & = (—w/U)* X 
(Ux/v,) of the order of 35 to 40, 144, = 4, P, = 1, and 
oT , a 

— = 0, the asymptotic profile is already approached 


OY 9 
for practical purposes. 
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APPROXIMATE CALCULATION OF HEAT TRANSFER 


A method of approximation for the heat transfer at 
the wall based on Lighthill!’ is easily formulated and 
First, define a scalar function \(y) 


* pu 
¥ -{ p dy (20) 
0 Pl 


which transforms Eqs. (33), (34), and (35) to 


dH _ l I pu dH 
dy P,dx . pi dd 


examined here. 
such that 


( pv) 


where 
H = c,T + (u?/2) 
If the assumption that the second term on the right 


side of Eq. (51) is small for P, number close to one, then 
this term can be neglected. Thus one obtains, 


(pv) (dH /dX) = (1/P,) (d/dX)|(upu/pi) (€H/dd)} (52) 
If Eqs.(5) and (33) are used, Eq.(52) can be integrated to 
(1 — T*)/N = (B/N) [1 — (u/U)]”" + 

1 — (u/U)? (53) 

which is the complementary solution of Eq. (41). The 

temperature gradient at the wall, from Eq. (53), is 
exactly that given by Eq. (48) with » = 1. 

If, in addition, the approximation that the velocity 


u is expanded in a series in \ and only the first term is 
retained of the form 


u/U = ar” (54) 


where a = V [(—pit1)29:]/[(pou0) U], Eq. (52) becomes 
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(pv) (dHT/dX) = [(upUa)/(P,p,)] (d/ddX) X 
[\ *(dH/dd)] (55 
which has a solution of 
(1 — T*)/N = (B/N)e2?/FPx”? + 1 — (u/U)? (56 


where F = [(up)Ua]/[P,pi(pv)] = constant. The 
temperature gradient according to Eq. (56) is 


[dT */dy|> = —BP,|[vo/(mC)] (11 To)? (57 
which is again exactly Eq. (48) form = 1. The two ap- 


proximations made for the asymptotic solution re- 
sulted in temperature gradients at the wall which are 
exact. Thus it appears that if heat transfer at the 
boundary for a more general problem is required, the 
approximations as given by Eqs. (52) and (54) would 
lead to values that should be close to the exact values 
obtained by more complicated methods. 


CONCLUSIONS 


An exact solution of the compressible laminar 
boundary layer over a flat plate with homogeneous 
suction or injection has been given for P, = 1 and an 
effective linear viscosity relation with temperature. 
The details for the suction case are worked out and are 
in excellent agreement with the results of the modified 
momentum method" except for the second derivatives 
of the velocity profiles which differ by about 20 per 
cent at the wall. A comparison of the first derivative 
at the wall is given in Fig. 6, which indicates good 
agreement in the suction case. This shows that the 
momentum method as given by references | and 9 
are sufficiently accurate for any practical purposes. 
However, any stability calculations should be based 
on the more exact values as calculated in this report 
and tabulated in reference 13. The characteristics 
of the injection case were calculated to a value of o = 
—0.2, and results are in good agreement with the 
momentum method of reference 2, although the results 
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LAMINAR COMPRESSIBLE BOUNDARY LAYER OVER 


of the suction case were in better agreement with the 
momentum method as indicated by a comparison of 
rigs. 5 and 8 for the first derivatives of the velocity. 
However, it is expected that the results of the momen- 
tum method differ considerably more from the exact 
calculations further downstream of the leading edge 
as separation of the flow is approached. Some further 
calculations are now being made, and results indicate 
a possible separation at 0 = —1.1. These results 


will be presented in a forthcoming report. 


The asymptotic solutions for the compressilbe bound- 
ary layer over a flat plate with homogeneous suction 
are completely derived for arbitrary heat transfer at 
the wall, arbitrary Prandtl Number, and a power law 
for the viscosity-temperature relation. Since the 
asymptotic solution can be easily obtained in closed 
form or by a simple quadrature, they are especially 
useful in determining the effect of the various param- 
eters. Closed form solutions are given for a linear 
viscosity-temperature law, and these results show 
that the effect of Prandtl Number on the velocity 
profiles is negligible and that the effect on the heat- 
transfer characteristics is large. 

The agreement of the exact solutions for suction in 
the first part of this report at large distances along the 
plate with the asymptotic solutions is excellent. It is 
noted that the asymptotic solutions are also valid for 
large values of the suction velocity at the wall, even 
with pressure gradient, compatible with the boundary- 
layer approximations. Therefore, these results can 
be used as first approximations in the practical calcu- 
lation of the characteristics of such flows. 
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Spanwise Variation in the Pitching-Moment 
Coefficient and the Center of Pressure Due to 
Various Basic Twist Distributions on 
Triangular Wings Having Supersonic 
Leading and Trailing Edges 


JULIAN H. KAINER* 


Ryan Aeronautical Company 
. pan. 


SUMMARY 


The linearized supersonic flow theory has been applied to de- 
termine the spanwise variation in the pitching-moment coeffi- 
cients and the centers of pressure due to various basic twist dis- 
tributions on triangular wings having supersonic leading and 
trailing edges 

The effects of the flow nonuniformity on the section properties 
produced by body upwash! and by geometric twist may be ob 
tained directly from the span loadings previously presented? and 
the results herein since the spanwise variation in incidence is 
known everywhere. To account for the effects of aeroelasticity 
where the spanwise variation in incidence is unknown, the span 
loadings? and the present results are applied to the superposition 
procedure reported on by Brown, Holtby, and Martin 

Design charts and generalized equations in closed form 
are presented for the section pitching-moment coefficients and 
the section centers of pressure corresponding to constant, linear, 
parabolic, and cubic twist distributions. Graphical procedures 
are indicated to determine the section bending-moment coeffi 
cients, the wing bending-moment coefficients, the wing pitching 
moment coefficients, and the corresponding wing centers of 
pressure due to the aerodynamic loads caused by the variation 
of angle of attack 

It was found that the effect of washout on the section char 
acteristics decreased with increasing degree of washout and 
with increased angle of sweep, and increased with Mach Number 
and the amount of washout at a given angle of attack. For ex 
ample, at the root chord of a 45° delta wing at 1J = 1.6 witha 
5° washout at 10° angle of attack the following results were 
noted: (1) a forward shift of about 1'/. per cent in the center 
of pressure, (2) a decrease of about 7 per cent in the section 
lift-curve slope, and (3) a decrease of about 11 per cent in the 
section pitching-moment slope over the wing with zero washout 
The net effect of this twist distribution resulted in: (1) a de 
crease of about 15 per cent in the wing lift-curve slope, (2) a 
decrease of about 19 per cent in the wing pitching-moment 
slope, (3) a decrease of about 25 per cent in the wing bending 
moment slope, and, (4) an inboard and forward shift of about 
2 per cent in the lateral and longitudinal centers of pressure, 


respect ively 


SYMBOLS 


Free-Stream Conditions 


V = velocity 

M = Mach Number 

B =VM-1 

m = Mach angle, aresin (1/.1/ 

p = mass density of air 

q = dynamic pressure, (1/2) pT? 


Received July 20, 1954 
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Wing Geometry 


5/2. = half-span 
Co = root chord 
Ce = average chord 
= half apex angle 
6 = angle between trailing edge and plane of symmetry 
m = 8 tan « 
Mo = Btand 
S = wing area 
x,y = Cartesian coordinate of system of axes with origin 


leading edge of root chord 


Analysis Parameters 


o = velocity potential 

¢, = horizontal perturbation velocity 

A; = are cos (1 — mt)/(m — 2) 

A» = arc cos (1 + mt)/(m + b) 

x" = x/cy evaluated at the trailing edge 

X,.». = section center of pressure 

Xe } ‘-2 x p c 

¥..», = Wing longitudinal center of pressure 

fg Bi Jes 

Vey = wing lateral center of pressure 

» ee = [(Yep) (b 2) ] 

a, = coefficients of power series expansion for the angle of 
attack 

c. = section normal-force coefficient per radian 

e. = section pitching-moment coefficient per radian 

Cy», = section bending-moment coefficient per radian 

a = wing lift coefficient per radian 

Cy = wing pitching-moment coefficient per radian (based 
on mean aerodynamic chord) 

Cn, = wing bending-moment coefficient per radian 

At) = conical functions from which the potential and normal 
force coefficients are determined 

Bt) = conical function from which the horizontal perturba 
tion velocities are determined 

et) = conical functions from which the section centers of 
pressure are determined 

é,(¢4) = conical functions obtained for JJ = 

Y = an arbitrary characteristic dimension 

Subscript 7 = index representing degree of twist 
Superscript * = represents multiplication by 8 


INTRODUCTION 


hes LOADING AND STABILITY of many classes of rigid 
wings in uniform supersonic streams have beet 
widely investigated both theoretically and exper 
mentally. The increased interest in high-speed pet 
formance has created a demand for more reliable design 
information; hence, a proportionate amount of re- 
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search has been directed toward some means of corre- 
lating available data with linearized theory and, par- 
ticularly, improving the present theory to account for 
such important effects as wing-body interference,‘ the 
consequences of boundary-layer propagation, and aero- 
elasticity.” > 

Empirical correlation* by the author of the span- 
load distributions on delta wings based on linearized 
theory with some supersonic wind-tunnel data has been 
sufficiently encouraging to warrant further theoretical 
study to improve the calculated distributions. 


* This classified information is available to properly cleared 


persons on request 


As pointed out by Beskin,! the presence of the body 
creates an upwash which varies along the span (and, 
in practice, along the chord of the wing It seems 
plausible that the effects of this nonuniform stream on 
the aerodynamic loading and stability may be deter 
mined upon consideration of a uniform stream in which 
is placed a wing warped to fit the stream lines of the 
nonuniform stream. The solution to this problem is, 
uniquely, that corresponding to a twisted and cambered 
wing in a uniform flow. Furthermore, the effects of 
aeroelasticity on the structural requirements may be 
determined by applying such an analysis to the super 
position procedure described by Brown, Holtby and 


Martin.® The noniterative procedure described therein 
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requires knowledge of the wing structural matrix and 
the aerodynamic loading properties? of the wing for 
various known basic angle of attack distributions. 

A systematic study” * ® is underway to provide the 
necessary aerodynamic data with which the solution® 
may proceed. The loading distribution for triangular 
wings having supersonic leading and trailing edges due 
to various basic twist distributions has been previously 
presented.” The purpose of the current paper is to 
provide a means of proportioning these loadings to the 
control points of the structural matrix. The spanwise 
distributions of the centers of pressure and the pitching- 


moment coefficients are presented here corresponding 


to the load distributions? due to known basic twist 
distributions for triangular wings having supersonic 
leading and trailing edges. From these results it is 
possible to determine graphically the section and wing 
bending-moments, wing lift and pitching-moments, 
and the wing longitudinal and lateral centers of pres- 
sure. Derivations have been initiated for the total 
lift, moments, and centers of pressure corresponding 
to the configurations considered here. Calculations 
are underway for similar configurations having various 
basic camber distributions and will subsequently be 
extended to include the class of triangular wings having 


subsonic leading edges and supersonic trailing edges. 
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ANALYSIS 


General Background 
The velocity potentials, the horizontal perturbation 
velocities, and the span-load distributions due to vart- 
ous basic twist distributions on triangular wings having 
supersonic leading and trailing edges were derived by 
means of linear theory.” The results may be concisely 
written, respectively, as 
(x/V) [(B'y)) x'*"] (@/a); = Aj(t) (1) 
(r/V) (B’x’)' (b,/a); = Bit) 2 


wey a C1 er yt* 
(C,/a); = At) (3) 


(3/8) [(mmo)/ (mo 


where A,(f) and B,(t) are the conical functions pre- 
sented in the Appendix of the original report.” 


Section Pitching-Moment Coefficient 
By definition the section pitching-moment coeflicient 
about the leading edge of the root chord ts 


*X trailing edge Ap 


Cc’ = A dx (da) 
“ *leading edge q 
= Min iGal (4b) 
Hence 
CCE = Dheg CHC Cee) (5) 


Substituting Eq. (3) into Eqs. (4) and (5) 


(x 16) [(mmo) (me — m)]i (1 x’)? +! X 
(Ge a) ;(¢ Gy.)* = AOa) 6) 


where 


such that 


0 <.<(f) S 10 


From Eq. (6) the resultant section pitching-moment 
coeflicient parameter is, by superposition, 


Cet) = 2 (Cul): (6/Cer.)"A (Sa) 
1 0 
= 2x' >" (C,,/a) ;aje,(t) (Sb) 
1 0 


Section Center of Pressure 
By definition: 
| section pitching-moment 


= : ae (Ya 
Co section normal torce 


, | ‘trailing edge Ap ‘trailing edge Ap 
= x ax ax 
Co x 


leading edge q “ *leading edge q 


ra , 
Xep 


Koy 


(9b) 


‘TE 
XrreOrre — f dx 
— (9c) 


Co Orr 
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(lower limit need not be considered since potenti] 
vanishes at a supersonic leading edge). 

Substituting Eq. (3) into Eq. (9c): 


OX TE 
| x!1+'4 (t) dx 
‘ 10a 


Xrp’t' A(t) 
e;(t 1Ob 


The results of Eq. (10b) are presented in the Appendix 
for values of 7 = 0, 1, 2, 3. 
From Eqs. (5) and (8), the resultant center of pressur 


is 


Wing Center of Pressure 

Derivations are currently being completed for the 
total wing lift coefficients, pitching-moment coeflicients, 
bending-moment coefficients, and the centers of pres 
sure corresponding to the various basic twist distribu 
tions. Until these equations are published the centers 
of pressure may be graphically determined as follows 


E-9. = Cec ay’ C,cdy’ 12 


Upon substitution of Eqs. (3) and (6) inte Eq. (12 


| (x’)? +44 (te (t) dy’ 


C7 ' 7) 


but (x’y’) may be related to the conical ray, /, at the 
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Fic. 9. Theoretical span-load distributions for a delta wing with 
and without twist « = 45°, AJ = 1.6, a,/ay = —0.5, ao = 10 
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Fic. 10. Theoretical pitching-moment coefficient distributions 
fora delta wing with and without twist « = 45°, JJ = 1.6, a;/ao = 
0.5 10 
trailing edge by” 
v’ = mo (mM — 2b) (14) 
y’ (my) — m)t/m(mo — ft) (15) 
and Eq. (13) becomes, finally, 


m | (my — t)~4-' A,(te,(t) dt 
eda v0 (16) 


(my — t)~*-' A,(t) dt 


Ina similar fashion the lateral position of the center of 
pressure may be determined for bending-moment in- 


vestigations by 


1 
e ( 


a ee (17) 


bite 1 
[ C,c dy’ 
e 0 


| t(my — t)~4-'A,(t) dt 


which becomes 


5 [my m) — 1] (1S) 


| (my — t)~*-* A,(t) dt 


Unswept Wings 

The results presented herein are also applicable within 
the root chord Mach cone of swept wings, hence the 
limiting case of an unswept wing may be determined 
by considering m = « These limiting cases are 
presented in the Appendix and are denoted by the 


barred quantity ¢,(¢). 
DISCUSSION 


The results of the analysis for the section pitching- 
moment coefficients and the section centers of pressure 
iS presented in the form of generalized equations in 


PITCHING-MOMENT COEFFICIENT 603 


closed form in addition to the design charts which are 
presented in parametric form and a numerical example 
for illustrative purposes. 

The section center of pressure parameters, ¢,(/) are 
presented in Figs. 1 to 4 for basic degrees of twist corre 
sponding toz = 0, 1, 2, 3, respectively, and the section 
pitching-moment coefficient parameters, ( 16 WA (bed), 
are presented in a similar fashion in Figs. 5 to 8. 
The results are plotted for several values of the leading 
edge Mach Number parameter, @ tan e¢, against the 
conical ray, ¢, evaluated at the trailing edge. The 
solutions are valid for all supersonic values of the 
trailing edge Mach Number parameter, § tan 6. 

The numerical results are illustrated in Figs. 9 to 12 
for a 45° delta wing at J/ 1.6 with a linear washout 
of 5° at an angle of attack of 10 Table I is presented 
to illustrate the manner in which the published results? 
are applied with the present results to obtain the section 
loading characteristics against semispan position. 
The span-load coefficient parameter corresponding to 
the twisted wing is compared to an untwisted wing of 
the same plan form in Fig. 9. At the root chord a 
decrease of about 7 per cent in the section lift slope 
multiplied by 6 is noted. The net effect of this twist 
ratio reduces the overall lift curve slope by about 15 per 
cent and produces an inboard shift of about 2 per cent 
in the lateral center of pressure. The spanwise vari 
ation in the pitching-moment coefficient parameter is 
similarly compared in Fig. 10. At the root chord a 
decrease of about 11 per cent is noted. The resultant 
decrease in the wing pitching-moment slope parameter 
is about 19 per cent with about a 2 per cent forward 
shift in the longitudinal center of pressure. These 
overall effects are predicted by some preliminary un- 
published* results for triangular wings at supersonic 


* To be presented when numerical results are complete 
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Fic. 11. Theoretical spanwise variation in the longitudinal 
center of pressure for a delta wing with and without twist « = 
45°, MW = 1.6, ag = 10°, a;/ag = 0, —0.5, —1.0 
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speeds which for a delta wing reduce to 


C.* a = | +- (4/3) (a, do) + (2/3) (d2/ady) + 
(2/5) (a3/do) +... (19) 
(Beg. /X le = 2/3; 3/4, 4/5; :5/6,.... (20) 

Cy*/ao = 4 + (3/2) (a,/ao) + (4/5) (42/ao) + 
(1/2) (a3/do) +... (21) 


The comparisons for the effects of twist on the longt- 
tudinal centers of pressure are presented in Fig. 11 for 
the linear twist ratios a; a) = 0, —0.5, —1.0. 

The spanwise variation in the bending-moment co- 
efficient parameter for a twisted delta wing is com- 
pared to the results for the corresponding untwisted 
wing. The effect of this washout is to reduce the root 
bending-moment coefficient parameter by about 25 
per cent. 

CONCLUSIONS 

The results of this analysis indicated that the effects 
of twist—1.e., body upwash— produced wing forces and 
moments which were greatly reduced in magnitude 
from those calculated for a wing having a constant 
angle of attack equal to the geometric angle of attack. 
Furthermore, small amounts of built-in twist may be 
used to reduce root chord bending-moments by a con- 
siderable magnitude. In general the effect of washout 
on the section characteristics decreased with increasing 
degree of washout and leading edge sweep, and _ in- 
creased with Mach Number and with the amount of 
washout at a given angle of attack. The results for a 
45° delta wing at JJ = 1.6 with a 5° washout at 10 
angle of attack indicated the following: (1) a forward 
shift of about 1!/. per cent in the root chord center of 
pressure and an overall forward shift of about 2 per cent 
in the longitudinal center of pressure, (2) a decrease of 
about 7 per cent in the root chord lift-curve slope and a 
net loss of about 15 per cent in the wing lift curve slope, 
(3) a decrease of about 11 per cent in the root chord 
pitching-moment slope with a net loss of about 19 per 
cent in the corresponding wing value, and (+) a decrease 
of about 25 per cent in the root chord bending-moment 
slope with an inboard shift of about 2 per cent in lateral 
center of pressure. 
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APPENDIX 


General Triangular Wings 
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Similar Solutions of Compressible Boundarvy- 
Laver Equations 


TING-YI LIf ann HENRY T. NAGAMATSU? 
California Institute of Technology 


ABSTRACT 


In this investigation, it is disclosed that the partial differential 
quations for the compressible laminar boundary layer can be 
transformed into a system of ordinary differential equations by a 
jirect procedure. This new procedure yields explicitly a differ 
ntial equation for the free-stream velocity distribution, 1:, 
On integrating this differential equa 


which 


llows similar solutions 
tion, the special forms of «, are found, and these results are corre 
ited with Stewartson’s results. As an example of the similar 
solutions, the hypersonic viscous flow over a flat plate is dis 
ussed. From the analysis of the numerical solutions, obtained 
m the Reeves Electronic Analogue Computer, the characteristics 
of the laminar boundary-layer flow with pressure gradient and 
heat transfer are discussed. It is shown (1) that on a heated 
surface, the local velocity inside the boundary layer can exceed 
the free-stream value when a certain favorable pressure gradient 
exists, (2) that the variation of Cy,/(1/2)Cy, for various surface 
temperatures with different free-stream pressure gradients can 
be accurately predicted, and (3) that the effects of heat transfer 
m the laminar separation phenomenon can be quantitatively 


studied 


SYMBOLS 


= coordinate axes 
= velocity components in the x, y direc 
tions, respectively 
= density 
= pressure 
= viscosity coefficient 
= coefficient of heat conductivity 
= gas constant 
= “similar’’ variable 
Prandtl Number 


total energy 


= temperature 

= stream function 

defined in Eqs. (20), (21), and (22), re 
spectively 

= Mach Number 

= ratio of specific heats 

= specific heat at constant pressure 

= speed of sound 

defined in Eq. (45 


the constant parameter in Eqs. (30) and 
(30a), see also Fig. 7 

C) = local heat-transfer coefficient 

Cy = local skin-friction coefficient 


Subscripts 0, 1, and w refer, respectively, to the stagnation 
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conditions, the free-stream conditions, and the surface condi 


tions 


I) INTRODUCTION 


6 ies SEARCH FOR “‘SIMILAK" SOLUTIONS of the com 
pressible boundary-layer equations can be dated 
from 1946 when Illingworth! examined the question of 
what forms of the variation of u,; with x would yield 
solutions with similar velocity and temperature distri 
butions for different values of x. His investigation 
demonstrated that such solutions exist only if the ex 
ternal velocity “; is constant, and then, as in incom 
pressible flow, the similarity variable is of the form 
y Vx. Later, Howarth,’ assuming that o = | and 
uw = (uo/7Ty)T, developed a theory of compressible 
laminar boundary-layer flow in which he again came 
to lllingworth’s early conclusion that no ‘‘similar’’ solu 
tions other than the one corresponding to “u; = con 
stant exist. In 1949, Stewartson* Mlingworth* 
independently demonstrated that the boundary-layer 
equations arising from a given main-stream velocity 
wu, of a compressible fluid could, by a transformation 
of coordinates, be converted into the equation for an 
incompressible fluid with a different mainstream ve 


and 


locity, 1). The relation between J; and y, is (1: 1 
thus, solutions known to yield “‘similar’’ velocity pro 
files in incompressible flow readily lead to solutions in 
compressible flow which are designated by Stewartson 
as ‘‘similar’’ solutions. The present authors, following 
Stewartson’s approach, have discussed the ordinary dif- 
ferential equations pertaining to these ‘‘similar’’ solu- 
tions; a preliminary report was published recently.' 
Strictly speaking, the ‘‘similar’’ solutions discussed 
in references 3 and 5 are solutions with similar velocity 
(w) and total energy (/7) distributions for different val 
They are, therefore, not similar as far as 
In the case x 


ues of x. 
temperature distribution is concerned. 
= constant, Illingworth! has shown that the momentum 
equation, the energy equation, and the continuity equa 
tion can all be reduced to ordinary differential equa- 
tions involving three unknowns, u, v, and 7’, which are 
functions of the variable yy In 
Stewartson’s approach, the variable v is eliminated by 
the introduction of the stream function y, which satis 
fies the continuity equation identically, and then the 


independent 


momentum equation and the energy equation are shown 
to be reducible to ordinary differential equations. 
Therefore, the introduction of y and // and the assump- 


tion of o = 1 and uw = (u/7y)7 are largely responsible 
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for the success of Stewartson’s theory.* In the pres- 
ent paper, it will be shown that, admitting these fea- 
tures into Illingworth’s initial work, the ‘‘similar’’ 
solutions can be deduced by quite a direct procedure 
without the need of introducing |, and transforming 
This new procedure yields explicitly a 
“similar” 


coordinates. 
differential equation for #, which allows 
solutions. On integrating this differential equation, 
the special forms of “, can be found, and these results 
are correlated with Stewartson’s results. Also, some 
new “‘similar’’ solutions will be presented. As an ex- 
ample of the ‘“‘similar’’ solutions, the hypersonic vis- 
cous flow over a flat plate is briefly discussed. Also, 
from the analysis of a number of numerical solutions, 
the behavior of a laminar boundary-layer flow with 


pressure gradient and heat transfer is discussed. 


(II) COMPRESSIBLE BOUNDARY-LAYER EQUATIONS 


The steady, two-dimensional flow of a compressible 
viscous fluid along a fixed boundary will be considered. 
With the coordinate axes of x and y along and at right 
angles, respectively, to the wall, the usual boundary- 
layer approximation yields the following equations: 


( Ou 4 ~ Op 1 a) ( | :' 
1 v = - (1) 
P\ Max * ” Oy dx | dy Yay 


0 = Op/oy (2) 
fe) re) Op 
ae— ¢. + ¢ cr ) =u-—+ 
o( ox * oy Ox 
(*) 4 3) (: ~) 
u , 
oy Ov \ Oy 
The equation of continuity is 
(0 /Ox) (pu) + (O0/Oy) (pv) = 0 (4) 
It is assumed that the fluid is a perfect gas so that the 
following relation holds: 
p = RpT (5) 


In the main stream, where the effects of viscosity and 

heat conduction are neglected, the pressure, density, 

temperature, and velocity are functions of x only. 
They satisfy the following equations: 

pil (diy ax) = —(dp, dx) (la) 

d a dp, Er i 1 dp, 


oy = came = ja) 
aor! te ide aa = 


fi = Rp,T\ (5a) 
The subscript ‘‘1’’ refers to the main-stream quantities. 
Eqs. (la) and (3a) imply that 
C,T, + (1/2)m,? = constant = C,7) (6) 
and thus, 
a;> + [(y — 1)/2]u,? = constant = a? — (6a) 
where 7) and dp are, respectively, the stagnation tem- 


* Of course, Howarth’s theory also contained these features. 
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perature and the stagnation velocity of sound in th, 
free stream. 
, " ») > ace 4 ’ ‘ 7 ‘ _ 
By Eq. (2) the pressure is constant, for a given valu. 
of x, throughout the boundary layer; consequently, for 


all values of y, 
p(x, yw) T(x, vy) = pilx)Ti(x 


It is further assumed that 


and with the introduction of 
H = GT + (1/2)n? 10 
Eqs. (1), (3), and (4) can be rewritten as follows 


u OU v Ou 


fi” Tie” 
G dT, Mo 1 oO ( - Ou 
~ + R 7 
T; dx ba To Pi Oy =) a 
u OT v Off Mo 1 Of, OH 
T Ox T oy Ty pi Oy Oy 


2) (;) oO (*,) _ =: (3) , 
ox \T) * ay\7) ~~ RT, dx \T , 


The boundary conditions are 
H = H, = GT, (14 
H= HM, = GT (15 


y> @ u= ty 

(III) RepucrBiILity TO ORDINARY DIFFERENTIAL 

EQUATIONS 

When the main-stream velocity is uniform, the terms 
involving (1/7) (d7\ dx) drop out from the boundary- 
layer equations. Then Eqs. (11), (12), and (13) can be 
reduced to three ordinary differential equations involv- 
ing the three unknowns, u, v, and 7), which are functions 
of the independent variable y/+/x. The possibility 
of obtaining a reduction from partial to ordinary dif- 
ferential equations when the main-stream velocity 
is nonuniform will now be examined. Let 


Eq. (13) can, then, be rewritten as 

v/T = —(C, RT;) (dT, dx)y — (Op Ox li 
and Eqs. (11) and (12) become 
Oy Ou =u ( C, dT; , Ow) | 
Oy Ox dy \RT) dx ax] 

7 Ss dT, LR ve 1 oO (7 3") 1s 
7, dx To pi Ov Ov 

Oy O/T =Oll ( C, aT, v4 =4 7 
Oy Ox = Oy \RT, dx * | Ox] 


wo 1 O (: : oe) 
x — 1 (19) 
Ty pi Oy oy 





Supp 
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tions, 
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When 


The it 
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Supp yse that 


Y(x, vy) = N(x)AK(n (20 
u(x, vy) = u(x) F(n (21) 
F(x, y) = M,G(n) (22) 


where 7 is the similarity variable. By definition, 
u = T(dy, oy) (23) 
Therefore, it is implied that 
T(On/Oy) = (u,/N) (F/K’ (24 


where the prime denotes d dn. After simple reduc- 


tions, Eqs. (1S) and (19) can be expressed as follows: 


dU (& dT; _ _KF'F _ 
“s RT. de | Nde/' Kk ~ 


C oF... wm tf FTF _\ 
= — x — — — — F (18a 
T, dx To piN? K’ \K’ 


7 +4. 
(¢ aw. = i 
- : T = KG = 
RT, dx N dx 


but 
C,T = H — (1/2)u? = MG — (1/2)u,2F? (25) 


therefore, Eq. (1Sa) can be rewritten as 
1 du Pe ( Ge ais l =| 
uy dx RT, dx N dx 
KFF’ _ 1 dT ( pr IT, G 
x - T; dx Fs Mu," : 


o uy F F ; 
R Le — ( a = r’) (1Sb) 
To pi? \K K 


By rearranging the terms in Eq. (1Sb), it can be reduced 


y-1 1 di =a 
re u, dx 
( C.. dis l I K FF’ 
aes 
RT, dx Ndx/} K’' 


0 ; FF ‘ 
R - : = : ( : r’) (18e 
T) piN*(du,/dx) K’ \K’ 


In Eq. (18c), let 


C, aT _1 dN _ ™m duy (: 7= 5 7 . Bye 
RT, dx © Ndx ~~ ™ dx J *) — 


Whence, it is found that 


1 dN m du; (m+ 2)y —m1 da, a 
= = — (2d) 
N dx u; dx y—-—1 a, dx 
The integration of Eq. (27) yields 
aera (98) 


N = Nou," (a) 


Where V, is an integration constant. Therefore, Eq. 


(1Sc) becomes an ordinary differential equation pro- 


vided that 


: ma (mt2 de 
Uy" (a, Pr 
dx 
= aod 
l | m\{1- M, 29 
) 
I 2 
where it has been assumed that .V R(o/To). Evi- 


dently, this same condition is sufficient to reduce Eq. 
(19a) to an ordinary differential equation. The 
momentum equation and the energy equation are, 
therefore, reducible to ordinary differential equations; 


they are, by letting F = K’, the following equations 
KK’ +k = (1/m) (Kh? = G (30 
G" + KG’ = 0 (31 
The boundary conditions are 
= K = K’=0 G =G, = T,/1 32 
n> K'=1 Ge=1 (33 


Here it is assumed that 7 = 0 or according to y = 0 
or ©. That this is the case can be seen from the equa- 


tion for n as follows: 


(1V) SUFFICIENT CONDITION OF REDUCIBILITY TO 
ORDINARY DIFFERENTIAL EQUATIONS 


It has been shown that if the main stream variables 
uy(x), Pi(x), ai(x) satisfy Eq. (29), the reduction from 
partial to ordinary differential equations is possible. 
Therefore, Eq. (29) can be regarded as a sufficient 
condition for reducibility to ordinary differential 


equations. Eq. (29) can be rewritten as 


| 


m-+ m+2 
_ : da, 


(d,") 


y¥-1 ( 2 ) — ne 
= (do~ (oo 
mpo ¥- I 


This differential equation must be integrated to yield 
a,(x), which then can be substituted into Eq. (6a) to 
Such 1;(x) functions would yield solutions 


(ag? — a,")’ 


find 1;(x). 
with ‘‘similar’’ velocity (wu) and total energy (//) dis- 
tributions for different values of x. 

Let 


(Q,/a =F (36) 


then Eq. (35) can be rewritten as 


, 1 ” 2 na 
aia = dé 
= E : 


The integration of Eq. (37) yields 
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ee 7 i) 2)4 
aQ-g°°*¢ wt «= 
/J 0 


at eee 2 sis ; 
A - dy x (38) 
mp yo 


where A is the integration constant. The integral on 
the left-hand side is an incomplete Beta integral which 
can also be expressed as follows: 


m+1 m+2 


*e 
} (i — &)” g dé 
70 
. » m-+-1)5 
y— 1 |; ! F (” — (m+ l1)y¥ 
Lm — (m + 1)7 y—1 
3 m—1-—my 


where F is the hypergeometric series, convergent for 
< l and m > 1/2. 


Cc 
S 


A combination of Eqs. (88) and (39) vields® 
| : 


~ 1 /a,\ 2m a (mV? 
7 ( ' , x 
m — (m+ 1l)y J\a 


[Tis + m)] [T(s)] = 1, 2, B,<-4. 


To find a,(x), it is necessary to invert the series given 
in Eq. (40); the computations are, therefore, in general, 


1.0 t a Fae 


0.2 
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rather complicated. A question of more immediate 
interest is that of correlating the present results with 
Stewartson’s ‘‘similar’’ solutions. 


(V) CORRELATION WITH STEWARTSON’S THEORY 


If m = 3 2, Eq. (35) becomes simply 
| ; 


1—6 2 


mea 2 
ay (da, dx) = —|{( y — 1)/spi a } 
which can be integrated easily to yield 
ay 3 — OY AM \s _ 
= = a xr * 1) 
3Po ) £2 


Therefore, by Eq. (6a), it is found that 


The value of m = 3 2 corresponds to Stewartson’s 


) ) 


8B = 2/3 = (2m'/m’ + 1) or m’ = 1/2; for such a 


case, Stewartson’s theory gives [Eq. (4.7), reference 3 


? 54 3 Bx\2—-* } 
“° = ao" 1 ~ ( + > 7 ‘ | ia 
: eal f a J 


) 





By letting B = (4 3po)ao Eqs. (43) and (48a 


become identical. In reference 5, it was pointed out 
that, for 
tions must satisfy an integral equation of the following 


‘“similar’’ solutions, the main-stream condi 


kind [see Eq. (34), reference 5]: 


- 37-1 i Pas 
m= —( mi 
ay ys 
l d 1 “x fa, 3 l . 
F Vdx (4 
pil,” dx J0 \aQo 


This last relation can also be written as 














Tw, 
FIG. |- NUMERICAL SOLUTIONS OF K'AND G (EQS. 55 AND 56) FORB =0.286 AND VARIOUS To Ss 
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lediate 0 
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w = 
ieee 
O T 7 
my - fy 6 Tw . 
To 
=1.0 
0.6 
+] 
t 
= =O 
° =06 
0.4 
) ‘ U 
{2 Ket 
0.2 
(G-1)/( De -1) = Ao 
r Ay A, 
#3 . | 4 u] 3 4 5 6 
rtson’s 
such a 
. . Tw 
nee 3 FIG. 2- NUMERICAL SOLUTIONS OF K'AND G (EQS. 55 AND 56) FOR f = 0.4 AND VARIOUS —*'s 
1a 1 fay\” bs 1 d » C8 fay? : then Eq. (35) may be written as 
= a (27) dx (44a 
os sai lle Jo \do uy" * (du,/dZ (1 mpo)ay" ; 16 
ba It can be verified that for m’ = 1 2, Eq. (43) actually Integrating Eq. (46) gives 
‘d out satisfies this last integral equation. It is perhaps inter- ; \1/2m=1 
aul , é , 2m — ‘ oe sit 
condi esting to point out that the integral equation, Eq. u, = ( Qo” Z) + constant (47 
owmg | 44a), actually constitutes the condition of reducibility mpr 
to ordinary differential equations. Its position in For m = 1 2, Eq. (46) can be integrated to yield 
Stewartson’s theory corresponds somewhat to Eq. (35) ; 
j = . =f. 1Y [(9 - 7 g 
in the present theory. My (constant) exp [(2 po) do Z| (4 
If, now, a transformation of coordinates, Z = Z(x), is The Z-transformation in Eq. (45) is analogous to 
1 introduced so that Stewartson’s transformation [Eq. (2.21), reference 3]. 
The special forms of #, in Eqs. (47) and (48) are anal- 
> m+2 —(m+1 ° . 
‘ , ° 7 : gous to those previously found by the present authors. 
Z= | (a," y—1 dx (45) ve : as I _ fi ; ‘ oe 
Jo [See Eqs. (26) and (27) in reference 5. | 
ay 0 
eee wa 
' U 
K's — 
0.8 “1 
\ 
0.6 
m'= -0.066 
0.4 
H-H, 
Hw-H, 
0.2 + , 4 To ' 
m=O QS ~ 
=-0.05 as 
ey Pe | 
aan 0 —— 
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FIG. 3- NUMERICAL SOLUTIONS OF K'AND G (EQS. 300 AND 3!) FOR 72 AND VARIOUS(m) S 
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(VI) DERIVATION OF SIMPLE 4;(x)'S FOR SPECIAL 
VALUES OF y 
For m = —[(2y — 1)/(y — 1)], Eq. (35) can be re- 


written as follows: 
(1 ) = : : d (**) 
Ao” dx \ay" 
9 54 l 3 27 l 
(y — 1)* Y 2 + — 
7 ayo ( ): y-1T (49) 
(2y —_ l pr 7 = l 


which can be easily integrated to yield 


; 9 5y—-3 
l 2 Pe ay 2(5 1) 
= 4 a’? {1 —- — a = 
ay — 3Ly — 1 Ay] - 


| 2 
rr DA ao%~! (x + x0) (50) 
Ties yi) 


Therefore, the special form of (x) becomes 


—1 4 y—1 


» a Y ~ Y 
— (22 : pm)” 3 Ay 57 (x ry xp)3 5 (31) 
oy = 
That is, the main stream represents a decelerating flow, 
atx = —xX, u; > ~, while asx —~ ©, 4; > m = O. 
For triatomic gases, such as COn, y = 9/7, Eq. (51) 
becomes 


1] - 
1/12, —(7/12 1/12 ~1) 
b= (i Po ay“! (% + xX) (5la) 


In the case of diatomic gases, such as No, Os, or air, y = 
7 5, Eq. (51) becomes 
9 1/10 (1/2) 1/10) - 
og is 20 Po ay y (x + Xo) (51b) 


For monatomic gases, such as He, y = 5/3, Eq. (51) be- 


comes 


- ’ 1/8 — (3/8 —(1/8 - 
uy = [(7/16)po] ’” ao (x + x9) )  (51e) 
The values of m corresponding to y = 9/7, 7/5, 5/3 
are, respectively, m = —11/2, —9/2, —7/2. For the 
case of an insulated surface, G = 1, Eq. (30) reduces to 


KK" + K’" = @(K” — 1) 8 = 1/m (30a) 


This last equation has been studied by Hartree.’ The 
solution was found numerically for a range of values of 


8. When 6 = —.1988, K’’ vanishes at 7 = 0, so that 
Ou/Oy vanishes when y = 0 for all values of x. There- 
fore, the values of m = —9/2, —7/2 may not be used 


because separation could already have taken place; 
however, in the case of a noninsulated surface, Eqs. 
(30) and (31) must be solved together; then separation 
does not always occur for these m’s (cf. later discussions 
of the effects of heat transfer on separation). 


(VII) Hypersonic Viscous FLow ON A FLAT PLATE 


The problem of the ‘‘strong’’ interaction between the 
leading edge shock and the viscous layer over a flat 
plate has been discussed recently by several authors.” * 
In reference 7, the present authors first deduced, by the 
momentum integral method, that the pressure on an 
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insulated surface behaves like x '/”; this prediction js 
in qualitative agreement with recent experimental 
findings. On the edge of the boundary layer it can 
easily be deduced from Eqs. (la), (8a), and (5a) that 


)\* 


pip: * = constant 5? 


The fact that this particular relation holds along the 
edge of the boundary layer makes it plausible to treat 
the problem of hypersonic viscous flow on a flat plate as 
a particular example of the “‘similar’’ solutions. Thus, 
Eq. (44) can be rewritten as follows for 1/,° > 1: 


“CF Ca MeO 
pi 27 || pi/prdx \p 
in () 3 vm (5: 
xX Jo \Po 


Since 
[x/(p1/po)] (d/dx) (p1/po) = —(1/2) 
it is deduced that 
— | l mn 7 - 
m’ = . B= = =- 54 
¥+ 1 m ma +3 7 


Therefore, for this particular case, the ordinary differ- 


ential equations are 
K'' + KK" = [(y — 1)/y](K" —G)_ (55 
G” + KG’ = 0 (56 


These equations were first obtained by Lees.’ Eqs. 
(55) and (56) with the boundary conditions of Eqs. 
(32) and (33) have been integrated on a REAC ma- 
chine with y = 1.4 and y = 1.67 for a range of values of 
T.,/T) < 1.¢ These solutions are presented in Figs. 
land 2. It is perhaps interesting to note that Lees in 
reference 8 treated only the cases of 7),,/ 7) = land 7, + 
T) = 0; for the latter case, an approximate equation 


instead of Eq. (55) was solved. 


(VIII) LAMINAR BOUNDARY-LAYER FLOW WITH 
PRESSURE GRADIENT AND HEAT TRANSFER 


The “‘similar’’ solutions are interesting because they 
are “exact’’ solutions of the boundary-layer flow with 
pressure gradient and heat transfer. A study of the 
“similar” solutions can provide some understanding of 
the effects of pressure gradient and heat transfer. In 
Figs. 1 and 2, some effects of cooling the surface are 
For a fixed 8 > O0—i.e., a fixed favorable 


exhibited. 
pressure gradient—these figures show that the skin 
friction, Tx = by(Ou/Oy),~, is lower when the suriace 


is cooled. This is because both yu, ~ 7, and (Ou/OY)w 
~ K’’(0) are smaller on a cooled surface. Tostudy the 

* At the outer edge of the boundary layer, strictly speaking, one 
should use the equations of inviscid rotational flow instead of 
Eqs. (la), (3a), and (5a). Therefore, this condition is correct 
only as a first approximation. 

+ The authors wish to acknowledge their sincere thanks to Dr. 
C. Gazley, Jr., of the RAND Corporation for his kind help 


obtaining these solutions. 
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FIG.4- NUMERICAL SOLUTION OF K' (EQ. 55) FOR =0.4, % 740 


effects of surface heating, a numerical integration of the 
ordinary differential equations, Eqs. (30) and (81), 
with the proper boundary conditions, Eqs. (32) and (33), 
has been carried out for a range of values of 7,,/7)> 1 
and for 8 > QOorg<0.* A typical set of data is shown 
This figure is particularly interesting be- 


in Fig. 3. 


* These results were previously discussed in reference 5. The 
wuthors are indebted to A. Presson of the Jet Propulsion Labo- 
ratory for his close cooperation in obtaining these numerical solu- 
tionsona REAC 
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FIG.S~"CRITICAL" VALUES OF Tw. vs. THE PARAMETER m' 
° 


cause it shows the effects of pressure gradient on the 
boundary-layer flow over a heated surface (7), 7 

2). It is observed that the effect of a favorable pres 
sure gradient is to increase the value of (Ou Oy), ~ 
K’’(0), consequently, the skin friction increases with a 
favorable pressure gradient; the effect of an ad- 
verse pressure gradient is just the opposite. It is also 
interesting to mention that when the surface is ex 
cessively heated, the local velocity in the boundary 
layer under a favorable pressure gradient may become 
faster than the free stream—for instance, when 7,’ 7) 
= 4.0, 8 = 0.40, a velocity profile (wu, = A’) of the 
form as in Fig. 4 is actually obtained. In Fig. 5 a curve 
showing the ‘‘critical’’ values of 7), 7) and m’ is given; 


for values of (7), 7) — 1) and m’ inside the shaded re- 
gion under the curve, the velocity ratio (uu, = K’) 


will not have values exceeding unity. Analysis of the 
numerical results of the ‘‘similar’’ solutions also pro- 
vides values of C,,/(1/2)C;, for various surface tem- 
peratures with different free-stream pressure gradients 
(Fig. 6). These data are obtained by using the follow- 


ing formula: 
(2C,,)/Cy, = —[1 (Ge — 1)] [G'(0) K’’(0)] (57) 


where 


; k, (Ol / Ov) x E Mi(OU/ OV) w i 
C, = : : CC. = : (oS) 
pill (G,, — ]) l 2) pitt," 


are respectively the local coefficients of heat transfer 
and skin friction. As expected, 2C,,/C;, = 1 on a flat 
surface with zero pressure gradient (m’ = 0, Fig. 6). 
For a flat surface with a favorable pressure gradient, 
dp,/dx < 0, m’ > 0, it is found that 2C,,/C;, < 1; the 
curves in Fig. 6 become steeper as 7),/7 > increases. 
For a flat surface with an adverse pressure gradient, 
dp, dx > 0, m’ < 0, it is found that 2C,,/C;, > 1; the 
negative m’ branch of the curve given in Fig. 6 is com- 
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TABLE | puted 
Solutions of KK” + K'”’ = BK” G) With Boundary Conditions: A = A’ =0,G = 7T./7) for 7 
G" + KG’ =0 ‘ K’ =1 iG a4 for »— The r 
ase, 
pB 0. 286 0.286 0.286 40 0.40 0.40 
m’ = B/(2 — B) 0. 1667 0. 1667 0. 1667 25 0.25 0.25 ind t 
K”(0) 0.540 0.670 0.765 5til 0.741 0.855 idver 
G’(0) 
0.409 0.505 0.520 188 0.517 0. 52¢ It 3 
7 To) l 
Cr, /(A/2C4,) 0.907 0.754 0. 680 0.870 0.698 0.6 transl 
7 To 0 0.6 1.0 Of 0 cont 
G 1 G 1 G 1* G 1 G l G f l 
” = K’ K’ K’ K’ K’ 
T w/To) l T w/To) | T w/ To) l 1 ] 1 7 ualit 
0.0 0 1.0 0 1.0 0 1.0 0 1.000 0 1.000 6 1.00 —_ 
0.1 0.05 0.95 0.06 0.95 0.065 0.950 0.948 0.069 0.944 0.077 0.94 
0.2 0.105 0.90 0.125 0.90 0.150 0.895 0.899 0.138 0.892 0.157 0.890 | 
0.3 0.160 0.855 0.190 0.85 0.215 0.850 0.850 0.205 0.841 0.23 0.838 studie 
0.4 0.210 0.805 0.250 0.80 0. 280 0.795 0.801 0.2 0.790 0.320 0.788 
0.5 0.265 0.760 0.310 0.75 0.345 0.740 0.754 0.335 0.740 0.369 0.738 Thus, 
0.6 0.315 0.705 0.365 0.70 0.405 0.690 0 0.708 0.396 0.690 0. 432 0. 68 sonar 
0.7 0.370 0.660 0.65 0.465 0.640 0.378 0.659 0.452 0.640 0.495 0.633 ey 
0.8 0.420 0.615 0.60 0.515 0.595 0.428 0.612 0.506 0.590 0.556 0.583 2 
OY” 0.470 0.565 0.55 0.565 0.545 0.480 0.568 0.558 0.542 0.606 0.533 B. 
ee 0.515 0.520 0.570 0.505 0.610 0.500 0.525 0.522 0.605 0.496 0.649 0.490 writte 
1.2 0.600 0.435 0.660 0.420 0.695 0.410 0.612 0.437 0.689 0.408 0.729 0.401 
1.4 0.650 0.350 0.730 0.340 0.770 0.325 0.690 0.355 0.761 0.323 0.796 0.319 
1.6 0.750 0. 280 0.795 0.265 0.830 0.250 0.759 0.282 0.822 0.252 0.855 0.247 
1.8 0.810 0.220 0 845 0.205 0.875 0.190 0.817 0.219 0.872 0.191 0.894 0. 185 In dis 
2.0 0.860 0.160 0.890 0.150 0.910 0.145 0.863 0.163 0.909 0.140 0.932 0.138 tis in 
2.2 0.900 0.115 0.920 0.110 0.940 0.100 0.902 0.121 0.938 0.100 0.953 0.097 
2.4 0.930 0.080 0.945 0.075 0.960 0.070 0.932 0.088 0.959 0.068 0.969 0 065 
2.6 0.955 0.050 0.960 0.050 0.970 0.040 0.953 0.060 0.973 0.043 0.981 0.045 
2.8 0.970 0.030 0.970 0.030 0.980 0.025 0.969 0.040 0.983 0.027 0.989 0.030 
3.0 0.980 0.020 0.980 0.020 0.990 0.015 0.980 0.027 0.991 0.015 0.993 0.022 [his 1 
3.2 0.990 0.010 0.985 0.010 0.990 0.010 0.986 0.018 0.994 0.009 0.996 0.013 F 
3.4 0.995 0 0.990 0.005 0.995 0.005 0.991 0.011 0.998 0.003 0.997 0.009 hea 
3.6 0.995 0.995 0.005 0.995 0 0.993 0.006 1.00 0 0.997 0.004 ») 
3.8 1.00 1.000 0 1.000 0.996 0.004 0.999 0.002 : | | 
1.0 0.997 0.003 1.0 0.000 wou 
‘2 1.000 6 
* [nsulated surface, G 1 is the solution 
3 
K 
( 
om ey > ey 
Ce (1 1) K"(O) 
To 20; : er T 
C nw (?%y )w 
——— h es yore } 
258 , Puy H, (2-1) a - 2m’ ; 
° m+ a 
| Ce, s sala pe | 0 
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puted for the case of a heated surface (7, 7, = 2 AB << By = 0(1 Asn — 0, from Eq. (59), it is found 
[he rapid rise of the value of 2C,, C,, for m’ < 0, in this that 
case, seems to indicate that C,, rapidly approaches zero 2 ee 
ind the boundary layer tends to separate under an lim A’’ lim BoA” A 
idverse pressure gradient of a certain magnitude. n—0 +0) A 

It is rather difficult to discuss the effects of heat | BAG As( kh’ 1) ] , 
transfer On separation since there are not yet a suffi ~ i K | o 


ient number of solutions with varying 7), 7)’s for 


ph dx > 0. Nevertheless, it is interesting to discuss 


This can be interpreted as follows 


jualitatively the separation of a compressible boundary K’'(0 [A’’’(0) Jago 
layer as follows ByAG — ApB(K’ 
Insulated surface: Egy. (30a) is the equation lim kK 62 
: 7 rae -() 
studied by Hartree.” When 8 -0.1988, A’’(0) = 0. 
hus, for 8 0.1988 Q, the boundary layer will Hence, it is easy to see that 
separate {pcr +r 
ee ; (0) — [A’’(0) 
2) Noninsulated surface: Eq. (30) can be re lim lac—or 
written as \G—0 AG 
, ; err} e 8 AB K’ 
A (1/K) [B(A G ' hee 59 lim | 63 
+() A AG kK 
\G—0 


In discussing the effect of heat transfer on separation, 


tis interesting to consider the following 


G 


| 


T AC 1 


AG | 60 


where 


The boundary conditions are such that A kK’ 0 
at n 0; therefore, from the right hand side of Eq. 
(63), it can be seen that the limit exists only when 





[his is the case of a surface with a very small amount 
1,A’’(O 
By) < 0, it is expected that forG = 1 + 


vanishes when 
AG, K’’(0 


f heat transfer. If for G = 


AB - Bi 


That is, the condition in Eq. (¢ 


AG 


(64) 


4) is necessary if A’’(0 


would vanish when 8 By) + AB, where AG < 1, and changes continuously with G. Therefore, for a heated 
TABLE 2 
Solutions of AA’ K = B( A G With Boundary Conditions: A = A’ =0,G = T,/7 for = () 
Cy + KG ) K’' = i. G S for y- 
0 HO7 0.40 ) 0. 02 0.10 
g/(2 0 50 0.25 0 0.01 0 05 
2 0 2.0 2.0 » 0 2 0 2 0 
K”(0 1.7386 1.4226 1.1154 0.480 0.415 0.120 
Gio 
0.612 0.588 0.560 0.480 0 460 0 39 
2 0.3520 0.4133 0.5021 090 iy) 28 
G ] G l G l G 1 G G I 
AK K K kh kK kK 
1 1 I ] 7 ] 7 1 1 1 1 ] 
0 0 1.000 0 1.000 0 1.000 0 1.000 0 1.0 0 1 00 
0.1 0.161 0.941 0.134 0.943 0.107 0.946 0. 046 0.955 0.042 0 958 0 O15 0 oF 
0.2 0.306 0.880 0.260 0.884 0.210 0.890 0 094 0.908 0.085 0 911 0.033 0 927 
0.3 0.433 0.818 0.368 0. 825 0.302 0.833 0.140 0.861 0.128 0 865 0 O54 0 S87 
0.4 0.541 0.756 0.465 0.765 0.386 0.777 0.187 0.814 0.171 0 S18 0 O7F 0. 847 
0.5 0.634 0.696 0.551 0.707 0. 463 0.721 0. 235 0.767 0.215 0.773 0 100 0 807 
0.¢ 0.712 0.637 0.627 0.649 0 533 0 H6e 0. 280 0.719 0 258 0.727 0 127 0.768 
a i 0.777 0.579 0.693 0.593 0.597 0.613 0.327 0.673 0.302 0 681 0 155 0.728 
0.8 0 831 0.524 0.750 0.538 0. 654 0.560 0.371 0.628 0.345 0 637 0 184 0 690 
0.9 0.875 0.469 0.799 0.485 0.705 0.508 0.417 0.582 0.389 0.592 0.216 0.650 
0 0.910 0.417 0.840 0 435 0.751 0.459 0.461 0.538 0.431 0. 548 0.249 0 611 
 : 0.959 0.326 0.903 0.342 0.826 0.368 0.545 0.453 0.514 0 466 0.317 0 536 
1.4 0.985 0.246 0.945 0. 262 0. 882 0.286 6.624 0.374 0.593 0 386 0.390 0 463 
6 0.997 0.180 0.972 0.194 0.924 0.216 0.696 0 302 0.667 0.314 0.464 0.394 
1.8 1.000 0.127 0.987 0.139 0.952 0.158 0.760 0.237 0.73 0.249 0 538 0 328 
2.0 0.088 0.995 0.097 0.972 0.112 0.816 0.182 0.792 0.193 0.611 0 269 
2.2 0.058 0.998 0.064 0. 984 0.077 0.863 0.135 0.843 0.145 0.689 0.215 
2.4 0.037 1.000 0.041 0.992 0.050 0.901 0.097 0 BRS 0.107 0.743 0 168 
2.¢ 0.023 0.025 0.996 0.032 0. 930 0.068 0.918 0.075 0.80) 0 128 
2.8 0.014 0.014 0.999 0.019 0.953 0.046 0.944 0.052 0.849 0 095 
3.0 0.008 0.008 1. 000 0.010 0.969 0.029 0.963 0.034 0.890 0. 068 
5.2 0.005 0.004 0.005 0.980 0.018 0.976 0.022 0.922 0 048 
+4 0.003 0.001 0.002 0.988 0.010 0.986 0.013 0.947 0 032 
3.f 0.002 0.000 0.000 0.993 0.005 0.992 0.008 0.966 0 022 
3.8 0.001 0.996 0.002 0.995 0.004 0.979 0 014 
1.0 0.000 0.998 0.000 0.998 0.002 0.988 0 009 
1.2 1.000 0.999 0.001 0.993 0 006 
t.4 1.000 0.000 0.997 0 004 
1.6 0.999 0 003 
1.8 1.000 0 002 
7 0.000 
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surface, AG > 0, Bo X (AG) > O; that is, B = B + 
Ag > (); this means that the boundary layer on a 
heated surface will separate under a comparatively 
weaker adverse pressure gradient than in the case of an 
insulated surface. This prediction is in qualitative 
agreement with the results obtained on the REAC 
for the case m’ < 0, T,,/T) = 2 (see Fig. 6). Thus, the 
boundary layer on a heated surface may separate for a 
value of @ in the range 0 > 8 > —O.1988. The effect 
of cooling is opposite to that of heating—i.e., the 
boundary layer on a cooled surface can sustain a 
stronger adverse pressure gradient. From a physical 
standpoint, the effect of heating the surface is to de- 
crease the density of the gas layer immediately adja- 
cent to the surface; this seems to decrease proportion- 
ately the kinetic energy of such a gas layer so that it 
can sustain only a very weak adverse pressure gradient. 


(IX) CONCLUSIONS 


(1) A direct procedure has been used whereby the re- 
ducibility of the partial to ordinary differential equa 
tions can be shown without the need of introducing 
VY, and transforming coordinates.* ° The reducibility 
condition gives the special form of (x) for which 
‘similar’ solutions are permitted. Correlation of the 
present results with Stewartson’s result is discussed. 
For m = 3/2, u(x) as given in Eq. (43) can be shown 
to be identical to Stewartson’s answer. 

(2) The forms of the 1,(x)'s for ‘‘similar’’ solutions 
are, in general, very complicated. Nevertheless in 
some special cases, simple forms of 1;(x«)'s are obtained 
in the present paper. 

(3) The hypersonic viscous flow past a flat plate can 
be formulated as a special case of the ‘‘similar’’ solu- 
tions. The pertinent value of 3 is (y — 1) y. Nu- 
merical solutions for this case are given in Figs. | and 
2. The data are also tabulated in Table 1. 

(4) Figure 3 shows the distribution of u/u, and [7 / Fi, 
inside a laminar boundary layer on a heated surface 
(1,7) = 2) for various pressure gradients in the free 
stream. The data are also tabulated in Table 2. 
The variation of the boundary-layer thickness under 
various pressure gradients is clearly seen here. 

(5) On an excessively heated surface, the boundary 
layer under a favorable pressure gradient can be so 
greatly energized that the local velocity in the boundary 
layer may be faster than the free stream (see Fig. 4). 
The data are tabulated in Table 3. 

(6) The variation of 2C,,/C;, for various surface 
temperatures with different free-stream pressure gra- 
dients is shown in Fig. 6. 

(7) The effects of heat transfer on separation are 
qualitatively discussed. The effect of surface heating 
is to cause the boundary-layer flow to separate at a 
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TABLE 3 
Solution of KK” + K’” = BK” G) With Boundar 
Conditions 


G’ + KG’ =0 K = K’'=0,G =TI,/7 


for 7 = 
K' = lG#=l1, lor »— 

0 G0) 0.605 
3 } T/T, — 1 105 
[ 7 1.0 Ch, 0 
. (1/2)¢ ~ 
K"(0) = 1.580 

rr G _ l 

’ . (7/te — 3 
0.0 0 1.00 
0.1 0.140 0.960 
0.2 0.275 0.880 
0.3 0.400 0.820 
0.4 0.510 0.760 
0.5 0.600 0.700 
0.6 0.685 0.640 
0.7 0.755 0.580 
0.8 0.815 0.530 
0.9 0.860 0.480 
1.0 0.900 0.425 
L.2 0.960 0.330 
1.4 0.995 0.180 
1.6 1.010 0.130 
1.8 1.020 0.090 
2.0 1.020 0.060 
2.2 1.020 0.040 
2.4 1.010 0.025 
2.6 1.005 0.015 
2.8 1.00 0.010 
3.0 0.005 
2.3 0.005 
3.4 0 


value of 6 in the range 0 > 6 > —O.19SS; the effect of 
surface cooling is to delay the separation of boundary 
layer to values of 8 << —0.1988. 


(8S) In Fig. 7, the values of m’ [from the equation 


B = (2m’)/(m’ + 1)] are plotted against 8. 
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When these last roots exist, it is necessary to investi 
gate which initial displacements (a;) from equilibrium 





In the following pages, the quasi-steady nonl near lead to oscillations about a9 and under what conditions —e 
P ‘ ° ° ° . ° e -_ : co \ e1g 
longitudinal oscillation of an aircraft with such lift and the motion ends at another value of S. bie 
ae : ‘ - ’ an oe 3 , : \ ing 
moment characteristics is examined in some detail. Phe lift slope curve assumed in this study is: R & 
nn ] . Xx 
The moment curve Cy(a@) has extrema at angles of C : u 
; ere a eee lear 
attack a, and a», and is represented as a cubic in a. ee 0 = nd } Radi 
rr . . . . (ao — a) (av — r) r< < , (( Kadit 
lhe lift curve slope C,,(a@) is constant up to an angle of : OT/ NG a si oe (0 
: . “1: : = 0 i im gyI 
attack a,,; it then decreases linearly until it vanishes pai ial Pv, 
‘ e % e Lilt c 
at the stall point a, and remains zero thereafter. Typical normalized lift slope and moment curves are slo 
The equations for the longitudinal oscillation of this shown on Fig. 1. Static 
aircraft are derived and expressed in terms of suitable If an aircraft flies at constant velocity LU’ along a Dam 
dimensionless parameters. Two approximate methods course inclined to the horizontal by an angle y(t), its ps 
are worked out to solve the initial displacement prob- equations of motion are C 
lem: The exact solution of the undamped motion C. 
a * * ° . . 4 lo ») [2 . a yo ~ ’ M 
equation is used as the first step of an iteration procedure (Wg) Uy = (1/2)pU°AC, I’ cos 4 (4 Flich 
when the damping is small. The exact solution of the . 0 . ss : ; 
; ; I 15 , : : (W/g)k*(¥ + & = (1/2)pU*Ac) Cul(a) + Air Sj 
equation with critical damping and a special nonlinear- [c/(2U)] [(Cue + C ~4tC (8 
2 ies : [€/\oU ) Uq 7 Vag) @ Wa) ( 
ity serves as the basis for calculations when the damp- igs For b 
ing is appreciable. Numerical examples are given to where IV is the aircraft weight, k is the longitudinal be 0. 
illustrate both methods. radius of gyration, g is the acceleration of gravity, curve 
pis the air density, A is the wing area, Cy,, is the damp 
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= . - Te mean aerodynamic chord. u 
The moment curve postulated in this investigation a ; mt di . I. . 
: . . . The following dimensionless parameters are intro €, 
is a cubic with extrema at a, a». It has the form: ; ing dimension parameters are 1 
duced: 
In ex 
Cy = Bla® — (3/2) (ar + a)a®? + 3aj;aa— K ep # 
| | r = ¢tV (ag) / (3k aay) ipo 
: . en (( negle 
where B, A are constants that will now be evaluated. a= a/a 
re . . . errors 
The constant B is related to the static margin ac, = , ; ; 
‘ Te Then, after some manipulation and under the two as is she 
since by definition: ; ‘ ; : : : 
: sumptions that a sin y <1 and (Cyg + City) wing 1S pro cumb 
dCy /dC,(O) = (1 Cre.) (dC /da) (0) = portional to C,,, the equations of motion yield prob 
(3Bajaz) re = —g/c (2) ? °a meth 
- seen rn ; rn eee 
a” + uC(a)a’ + M(a@) = pe C(a) da 10 ind t 
The aircraft is trimmed at an angle of attack ay by are critic: 
application of a tail moment just sufficient to cancel where Bo 
the wing moment. If one examines a fixed control . " phase 
5 B 3 : M(a) = (&@ — &) a + 1@ — (3 2) (1 + @)/@ + : ; 
maneuver, that tail moment remains essentially con- (3/2) 1 4+ rie 1] &) 3 
ee Ps _—” = te , 4% QA) TT Aah) 1 
stant throughout the maneuver, and A gives its con- the te 


tribution to the total moment at any time; A is there re wes ; ee ; 
; . ‘ i le rhese coefficients include contributions from both wing 
fore selected to make Cy(@) vanish when a = a. The and tail. The wing contribution is assumed proportional to lift 
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QUASI-STEADY NONLINEAR 
a = 0 a a 
@)/(1 a a a l 
(delta aircraft (llb 
0 l<@ 
Fe aircraft with tail) 
Mi ( Va + 1] SR Gay 
Li | lle 
2k Cle al “a 
te 2k | (C MQ t. a 
1 c? (2k?) | (Cua + ¢ Cheol | 


[he important dimensionless numerical coefficients are 


therefore the damping coefficient » and the coefficient e, 


which measures the effect of flight path curvature. 
In general, at high speed and high altitude, u is small, 
Two typical aircraft 


iid eu” is small in most cases. 


re illustrated below, and yu, € are computed for them: 


(A B 
30,000 Ibs. 30,000 Ibs. 
200 sq.ft. 300 sq.ft. 


Weight 


Wing area 


R & Shape 2.31 (60° delta) 1.00 (45° swept) 
\ean chord 9.2 ft. [2.7 ft. 
Radius of 

3.9 ft. 6.2 ft. 


gyration 


Lift curve 
slope Cie 3.0 radians ~! 1.0 radians 


Static margin 1.8 ft. 0.4 ft. 
Damping 
coefficients 
co —0.75 -2.5 
Cye —0.60 —2.0 


Sea level 
100 knots 


35,000 ft. 
600 knots 


Flight altitude 


\ir speed 


For both of these aircraft, the stall angle is selected to 
be 0.40 radian and the first extremum of the moment 


curve is at 0.16 radian. It is then found that 


(ao), = 0.136 radian (ao)g = 0.0079 radian 
Ma = 0.107 ig = 1.04 
€4 = —(0.132 é€p = —0.074 


Inexample (A) where u, eu? < | the solution is expanded 


in powers of w and only the first term is retained; the 
neglected terms being of order yu? introduce negligible 
errors. While a scheme for finding higher order terms 
is shown, the calculations involved rapidly become so 
cumbersome that if the first term is not sufficient, the 
problem may be better handled by numerical or analog 
In example (B), the term in eu” is neglected, 
and the equation is solved by relating it to a special 


methods. 


critically damped nonlinear motion. 

Both solutions are conveniently examined in the 
phase plane whose coordinates are @ and ff. If @’ = 
'(@) is introduced into the fundamental Eq. (10) with 


the term eu? neglected, it becomes 


+ MW(a) = 0 (12a 


ff’ + wC(a) 


where 


PITCHING OSCILLATIONS 619 


r= > 12b 


with the boundary condition /(& Q. The point &, 
represents an initial disturbance from trimmed condi 
tions of flight, which may be due either to entry into 
a sharp gust or to a displacement of the control sur 
faces. The problem is then to trace the history of the 
motion as a function of trimmed angle of attack & and 


initial displacement a 


3) SLIGHTLY DAMPED NONLINEAR OSCILLATIONS 


If uw is sufficiently small, Eq. (12a) can be solved 


formally by expanding /(a@) into a power series in y, 
substituting this series into the equation and equating 
coefficients of like powers of u. The first term of such 
a series is the solution of the undamped problem /y(a 
which is now discussed in some detail. It satisfies the 
equation : 


(lia 


"a a ; . 
T\ : hy( a 0 (13b 
x. f 


(a, a 2 M(@aa = b 
J & 
= —(1/2) (&€ — &@ ;a? + 
, 14 
a, — 2(&@ + 1)|a@° + 
|a;7 — 2a(a@ + 1) + 6&)Ja + 
a; — 2a + 1)a° + 6aa& — 1K (Go) |f | 


The direction of integration must be so chosen that 
the integral is negative. If at &,, 1/(&) is positive, the 
direction of integration is negative and vice versa 
This requirement is equivalent to the physical con 
dition that motion proceeds in the direction of the 
spring force. 

Since the moment curve ./(&) in general has two 
stable focal points S;, S;, and a saddle point S»,* the 
trajectory field is divided into three regions R;, Ro, Rs; 
by the double trajectory through the saddle point S» 
That loop is given by 


/,.° = -2f M(a@) da | 
7 S - 
‘ . 6 aan (15) 
(1,/2)8,2(8.2 + 4).S. - | 
[(a@: + 7) 2j} +608, — a) (S. 1)) 
where Bs = &@ — So, the distance from the saddle point 
S.. It follows from the geometry of 1/(a@) that a < 
So < 1 so that (So — &) (S: — 1) < 0, and f<? exists 
when 
(1 + @ — Ss) — 
Vil + & 2S; — 6(S a) (S 1) < 
a (l+a,-—S + 
V (1 + & — 282)? — 6(S2— &) (S2— 1) (16) 


In particular fs) vanishes at the two extreme points @,, 


&,, defined by (16). If a < & < S, the trajectory 


f(&;) is a closed loop in region R;, and the motion is a 
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periodic oscillation about the focus S;. If So < a; < At S; and S;, 7 takes the values predicted by a linear ” 
@,,, the trajectory is a periodic oscillation about $3 in perturbation of the moment curve about those points we 
R;. If either &; < &, or &; > &,,, the motion is a large For a given value of @ [aircraft (A)], 7(&;) is plotted ep 
amplitude oscillation about S: in R: which is exterior on Fig. 2. wh 
to R,, R;. The polynomial }6)(@&&@;) has four real Returning now to the differential Eq. (12) which be - 
roots— two each in R, R3, if inequality (16) holds for comes after one integration: 

a@;, and two real roots in Re» if (16) does not hold. - co" 

If 1/(&) has only one root, it is a focus and the motion f? = bo(&) — Qu | fC da (2 le 
is not essentially different from linear motion. That — 
case is not treated here; it has been discussed by one defines a sequence of functions f/,, as follows: 

Oswald.° - 

If the four roots of the equation ))(@;, ®) = O are de fn? = bo(&) — 2p | fe aa (29 It 
noted by r;(&;, &), one writes: wits r, 
fo? = —(1/2) [a? — (n + ra + nro] X If the sequence /, converges, it converges to the re ws 

la? — (ns + 18+ 1) (17) quired solution. Now, one has 
The negative sign is due to the fact that either 7, < f. = Whim) [ a. Joc, bol a ia 
@< r2in R; or rs < & < ryin Rif all ry; are real and r by (GJ a D 
to <3 <3; OF 1) < & < ro in Re with n, ro real and 73, V bo[l — (2ubi/by) | (234 D 
rs complex conjugates, so that in all cases the product ; : : 
L iemaeliee. and taking the first term of the binomial expansion 
. Substitution of Eq. (17) into Eq. (13a) gives the solu fy, = V by [1 — p(di/bo) ] 23h “ 
tion as follows: 
: a Similarly for fo, f;, after integrating by parts and ex n 
i. i) da panding up to ’, u* one finds: 
ai —(1/2) [& (ri + r2)& + rire] X 7a b "4 
\ [a> — (r3 + 74)& + rgry] (1S) fo = mn — 2u | Cx i(1 = u') ia Ds 
J a Ay 
where one of the roots r; 1s the initial value @; and &@ 7 , 2 /p.2 
varies between @; and its associated root. This integral = % ak — 4p a « (; os 
equation may be inverted in terms of Jacobian elliptic b by \2b ; 4a | 
functions to give: / a ia) etc 
J/J@ r | ¥ 
a(t) = £6 —b oe cn k —~ a T; = i _ 2 . ” 
0 ) \ 5 on = MW, bh , = 0 [1 (by bo) + w?(be bo) | 


de : { ) a 1 — Xo ae J, _ ac . pd, # (x 
\' . ( I \' \; * \ ig J & | by T bo \ 2D, 


= he 


Ay 
= 2 "# bb,’ ) 
mi b fi on fy) ie | Ov ga) by" | aay 
: ,eory J & yn 
Az — Na ES Wr — 1) \ \y f 
b bs 3 [bibs 
where 1 oO Me 102 
= Vd) 1 - er - < + (241 
™ | . by “ bo b (= 
Ay do = 2G » wt 
* bob’ * 
[(n1 + re) (r3 + 74) — 2(nre + r3r4) F j da 
J & Hy 
2V (rn — rs) (14 — Pe) (1 — 14) (rg — Pe) | 
(rf, — 4) ='S by [1 iia u(d, by) + u"(do by) = 
(19b) u*(b3, bo) | 
and generally 
a,b = V (nre — Aa, 237s) /(1 — Ad, 2) (19¢) ” ‘ 
; . —e : : s Ds 
and cn(or, k) is a Jacobian elliptic cosine function of f= Vb dX (—w) F 
, <e : " ro 
modulus k. If the roots are positive, it can be shown 
that \;, A» are both real negative numbers; if two roots - ; =< bib, —1 
. ; ' b, = da + oa 
are complex conjugates, Ai, A» are of opposite sign. Jai do bo 
The period of the motion is bob, _» Bey 
j + + € oy (20a) 
] Ny 20h 
To = [(Ae _ A) ia- b)] \ [2\1(\; —1))/1—- A2) P 4 
K(jA2/A1}) (20) where 
The period has three singularities at @,, &,, (the dis- = 2f M da, bh, = i) cri d&, & = : : preg 
~ - O } 


placement &; is chosen so that the trajectory is the 


singular one through S.), and S», the nodal point. (25b) 





’ a linear 
‘7 points 
S plotted 


vhich be 


ion 
(23b 


and ex- 


(24b 


n) 
') 


25b) 





QUASI-STEADY NONLINEAR 


This sequence converges nonuniformly if uw is sufficiently 
small since b,/by is of order n. Special care must be 
exercised in the region where by) < yw (near &, and &;;) 
where the expansion process breaks down, and one must 
iterate on Eq. (12) directly. 

The time history of the motion in terms of the initial 
conditions is calculated by carrying out the integration 


13b) with f(@) given by Eqs. (25a) and (25b). 


f da [ da 
T= = (26 
ai f(&) JaVby dX, [(—H)"bn/ bo] 


/a 


It is convenient to take the reciprocal of the series in 


Eq. (26) as follows: 


%a 1a 
T= > u"D,,(&) in (27a) 
J ai 0 V bo 
D l 
D, = b,/6 
| °& fF 2} , 
a Db,? Ja Do? 
by [ by,’ l Kew (27b) 
, = 1¢ la <«dD 
L Dbyp? Jai by? veil bo? J ai Oy ea 
bs? 2b," [ by"b,’ 
D, = a ca — 
; 2by* bot Jai do" 
l [ bsb; 
da 
bo? J ai Dg 


etc., and the integral (27a) is evaluated term by term. 
The first term is given by the preceding analysis of un- 
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damped motion and the other terms are corrections to 
it so that 


T= %mt+ } i br (27¢c 
n l 
a la 
— [ D. da (27d 
J &i V bo 


and since the parameters D, are rational algebraic ex 
pressions, the corrections 7, are expressible in terms of 
elliptic functions. 

The general structure of the trajectory field in the 
phase plane (Fig. 3) is best studied by considering the 
four critical semitrajectories through the saddle point 
So. The trajectories 7), 73, which start from S) at 7 = 
— « toward 5S; and S; with positive slope f’ = — MWC. 2 
+ V wC?/4 + A’, always lie between the @ axis and 
the undamped trajectories 7\)73) through S:. They 
must therefore converge to S; and S; or a closed loop 
surrounding 3. 

7;, which reach Ss at r —~> + 
V w?C?/4 + MW’ 
coming from outside, divide the phase plane into two 


The trajectories 7», 
with negative slope f’ = —yuC/2 
spiral-shaped noncommunicating domains R, and R;, 
which contain the stable focal points S; and S3, respec 
tively. Since the damping function C(&@) is everywhere 
nonnegative, the trajectory through any point of R, 
converges to S; in a finite period of time, and the tra 
jectory through any point of Rs; similarly converges 
to S; or to a closed loop including 5;. 

Because of their spiral-like nature, the trajectories 
7», Ty intersect the & axis at a doubly infinite sequence 


of points. 7», intersects the @ axis at points a 





12 














" fs Se | [Ss 
| | 
| | 
| | | 
} } 
| Pt | 
| 
| | 

















———EE 





























en a 





| 
; 
| 
| 
| 
| 
| 
| 





| 
6 8 
a; 


1.0 1.2 1.4 





Fic. 2. Variation of period with initial displacement for pitch of aircraft (A) 
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dy, Go, dz... to the left of S; and at points @ = by, be, 
hb, ... to the right of S;. 7, intersects the &@ axis at 
points & = ¢, é, &... to the left of S; and at points 
@ = di, do, d3... to the right of S3. 
are related by the inequalities: 


These sets of points 


a>i>a>&m>a>6...\ (28) 
d, ‘. b, a ds < bs <x ds < bs oi sl . 


Any trajectory with initial value @; defined by 


Ge > es > cy | 


. (29a) 
b, < Qj < dn+i { 


converges to $3; any trajectory with initial value @; de- 
fined by 


Gig oo My e Cy if (29b) 
4.<@:< 6, 4} 

converges to S;. The number of large amplitude os- 
cillations of such a trajectory is m. Note that a tra- 
jectory does not always converge to the nearest stable 
singular point. The case (29a) in particular, corre- 
sponds to the observed pitch-up phenomenon observed 
in transonic flight. 

In view of this discussion, it is interesting to estimate 
the location of the critical points @,, b,,, &», d, Which are 


defined by: 


f(An) =@ A, re); mee (30) 


From Eq. (21), one obtains, for ay: 


ey 

by(a,) — au | fCda=0 (31) 
Se 

If uw is sufficiently small, the integral in Eq. 


small and one may write approximately 


bo(G1) = bo(&) + (4, — &) X 
(db) /da@) = 2M(&@,) (a1 — &7) 


(32) 
%al 
so that up to the first power of u 
; _ m al ; 5 . ™ 
ai —- @; = meee i fo C dx (33) 
Similarly 
m ail ; 
— a1, = foC da (34) 
d QyI 2M (a7) f : a 


since a; is outside @, while d; is outside @;;.. The other 


points are defined by 


ay =, ») mas a foC d& 
h ied _ oC : OL 
a dy ee [s C dx 

(35) 
qi — @e = om foC da 
da — 1 = ony ey I “om 


ClO. 2.45 
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From Eqs. (34) and (35) the entire pattern of critica] 
points da,, b,, €, d», is determined approximately jn 
terms of S.(&), and the final outcome of any motion 
characterized by an initial trim @ and an initial dis 
placement @; is determined from inequalities (29a, 
29b). This is illustrated on Fig. 4 where the abscissa 
is @& and the ordinate is @;. Every point of the 
(@;, &) plane specifies a trimmed attitude and an 
initial displacement and therefore defines an oscillation 
The curves a, (a), 5, (a), ete., give the intersections of 
the trajectories 7°; with the @ axis for all values of a 
and, consequently, the ranges of @; for any value of 
& which lead to motions in R; and in R;. Therefore. 
by plotting the critical points a,, 5,, €,, d, as functions 
of & and cross hatching the regions where (29b) holds, 
one can predict the outcome of any motion. If the 
specified initial point (@;, &) is in a cross-hatched region, 
the motion converges to .S,; if it is in a blank region, the 
motion converges to S;. Furthermore, the number of 
stripes to be crossed in going from &; to the equilibrium 
point along a vertical line is the number of large-scale 
oscillations. 

A numerical example of slightly damped nonlinear 
pitch-up is shown on Fig. 5. The trimmed angle of 
attack is @ = 0.34 [aircraft (A)], and the initial dis- 
turbance is defined by &, = 0.15, which corresponds to a 
sudden normal acceleration of 0.55g. It is seen that 
although the nearest stable equilibrium point is S; = 
& = 0.54, a sufficient initial acceleration is imposed 
to cause pitch-up ard an oscillation about the stable 
point & = S; = 1.29. Since the oscillation eventually 
takes place completely in the stalled region (@ > 1) in 
which the damping vanishes, the final state is a half- 
limit cycle instead of a subsidence in the neighborhood 


of S3. 


(4) HEAvILy DAMPED NONLINEAR OSCILLATIONS 


> 


The method described in Section (3) converges in 
However, it converges more slowly 
and the labor re- 


some range uw < po. 
as the damping increases (u — po) 
quired to get a solution increases. It is therefore 
worth while to look for a different approximation which 
retains the characteristics of the nonlinear equation 
and is based on a known exact solution. 

For purposes of convenience the trimmed equilibrium 
point will now be taken as the origin. The governing 
differential equation is of the form 


a” + pa’ + M(a) = 0 (36a 
or more specifically 
a” + pa’ 7+aa>+ba+t+c=0 (36b 


where the coefficients a, b, c are given in terms of the 


airplane characteristics by 


b = 5a 
C = —&l[&? — (3/2) (1 + &) + 3a] 
Note that in this investigation, one assumes C(&) = |. 
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If C(a@) changes its value, the approximations to be 
presented can be corrected by replacing yr by 
u J C(adr. In this substitution it is most important 
that the period and phase of the approximation be as 
accurate as possible. 
If the dependent variable is changed to 8 where 


B=a-— %&% 
the governing equation becomes 


B” + pB’ + B® + (8a + a)B? + 
(3a + 2ax + b)B = 0) 


(36c) 


This equation can be reduced to the canonical form 


by the substitution 8 = ye “’~, which gives 


y” + [3a + 2aX + b = (u> 1) Jy + 


J oF 


(3a + a)y*e ”’""" 4- ye = 0 (37) 


with initial conditions 


yo = ¥(O) = & — XH 


Yo" = y’(0) = £'(0) + (yu 2) Vo = | (oS) 


a’(O) + (u/2) (&: — &) 


Equations (36) or (37) can be integrated exactly’ 
if the coefficient of y® is zero and if the damping co- 


efficient takes the value 


uw = (3/V 2) V 3a? + 2am + b (39) 


This value of » can be related to the critical damping 
ratio ¢, if ¢ is defined as ¢ = yw/2w, and w,” is the 
coefficient of 8. Making the substitutions it is seen 

‘ o Is /5 . 
that ¢ = 3/2V/2 = 1.06 
The exact solution is then 


#T on (Ay V 3&2 + 2am + bX 
+ 2ado+b oo ¢; 1 \ 2) 


B = Ave” 
(40) 


e-(r 2) V 3do? 


and the constants A» and ¢ can be calculated in terms 
of the initial conditions from the formulas 


Ay —_ Bo(1 a ) (Bo’ + Bo) [B0?(3. a" + 


2am + b)}}2)'* (41a) 
¢ = en=![(Bo/Ao), (1/7 2)] — 
Ay V 3&7 + 2aa +b (41b) 


If the equations were linear—.e., if the coefficient 
of y were extremely large in comparison to unity—Eq. 


(37) would have the solution 


y = BoV 1 + [(a0’/Bo) + (u/2)]?/w? X 


cos }wr — tan! [(ao’/Bo) + (u/2)]/w} (42) 
which can be written in terms of 6 as 
B = pe *?V 1 + [(ao’/Bo) + (u/2)]?/w? X 

cos }wr — tan~![(ao’/Bo) + (u/2)]/w} (43) 


where 


w®W = V 3a"? + 2am + b — (u?/4) 


AERONAUTICAL 


SCIENCES SEPTEMBER, 1955 

The linear solution [Eq. (43)] is valid for all y as 
long as w is real. The solution to the nonlinear ecya. 
tion [Eq. (40)] as noted before is valid for only one par. 
ticular value of w(u = 2,¢), which corresponds t 
heavily damped motion. 

Examination of Eqs. (40) and (43) indicates certain 
similarities. they both have exponential multipliers 
and have a phase angle introduced to allow for the 
effects of damping upon the motion at the initial point 
Equation (40) will be useful in verifying the accuracy 
of the approximations. . 

The simplest approximation for the solution of Eg 
(37) consists in replacing the time dependent exponen- 
tials by constants, so that if 


wT /2 


7. = © (44 
where 7’ is the period of the motion or its linear approxi- 
mation, then Eq. (37) becomes 


y” + [Ba&o? + 2am + b — (u?/4)]y + 
(3a + any? + n*yi? = 0 (45 


This is an undamped equation and its solution may be 
found by the method of Section (3). Combination of 
the solution of Eq. (45) with the exponential gives 
B =e #l*y (46 

This approximation reduces the effectiveness of the 
damping in the earlier part of the cycle and increases 
it in the later part. 

It was demonstrated earlier that the nature of the 
solution depends upon the roots of a polynomial, in 
this case the quartic 


bp = Ci — [Ba ? + 2am + b — (u?/4) ly? — 
(2/3) ¥1(3G + a)y®? — (17/2)y4 = 0 (47 


where 


Ci = wo? + (9717/2) y04 + (2/3) 1(3& + a)yo? + 
[3a + 2aa& + b — (u?/4) |p’ 


It is sufficient to say here that if the quartic has four 
real roots, vi > Yo > y3 > ys, the coefficient of 7 in the 
solution is *® 


[y/(2V2)] V (91 — ys) (v2 — ys) (48 


If the quartic has two real roots y; > y2 and a complex 
the roots of (y — r)? + S?;S>0 


pair the coefficient 


of 7 1s° 
(Y1 2) [s(y1 = FF ky 1 — k?] (49 


where k is the modulus defined in Section (3). 

Note that 7; appears above explicitly and will appear 
in the expression for the roots. For example, if the 
cubic term in Eq. (47) is zero, the solution can be 


written as 


where A, and ¢ are as previously defined and 
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STEADY 


The objection to the use of a constant y; can be re- 
moved by replacing 7: by the exponential e 
functions y(7) no longer satisfy the differential equa- 
tion, but the approximation to the solution is im- 
proved. The 
solution as w approaches zero and behaves properly 


approximation approaches the exact 
for large time when the exponential term is small. 
This approximation will be useful if the equation is 
nearly linear or if the motion decays so slowly that the 
oscillatory terms have a smaller period than the char- 
The 


could be either a moderately damped motion where the 


acteristic time of the exponential. latter case 
linear and nonlinear terms are of the same order or a 
strongly damped motion with a large nonlinear term. 
The exact solution, Eq. (40), the first approximation 
y1 = constant), Eq. (50), and the suggested approxi- 


mation 


V 3a? + 2am +b — (u2/4) + 0” Apter + ¢: RI 


have been calculated for a motion whose damping 
would be approximately critical if the spring force were 
linear. The results are compared on Figs. 6, 7, and 8. 
When the nonlinear terms are either small or large 
compared to the damping and linear spring terms, the 
relaxation of the restriction to constant y; improves 
When 


the nonlinear terms are of the same order as the damping 


the agreement, which is considered satisfactory. 


and linear spring terms, the agreement is less satisfac- 
tory. 

In the general case the roots of the quartic are com- 
plicated functions of time. It seems more practical 
to compute them numerically either by Newton's or 
Lin’s method at several values of time. The variation 
of the roots with time can be approximated as suggested 
in the Appendix. Once the roots are known, the re- 
quired constants can be computed as functions of time 
and the values of the elliptic functions found from a 
table.’ The amount of work seems to be independent 
of u, depending only on the time interval selected. 

As an application of the methods described above, 
consider airplane (B) of Section (1). In normalized 
form the equation of motion is, if the damping is ap- 
proximately constant, 


” 


& + 1.04a’ + a — 2.10a? + 
117a@ — 0.115 = 0 (53) 


with initial conditions 


—(0.061 54 
> (.>- } 


= —0.75% = 


ie =e { 


These conditions correspond roughly to a —0.75g 
The reduced equation about the equilibrium 
point @ = 0.12 is 


gust. 


v" + 0.438y — 1.74y2e~ °°?" 4+ yie~' 7 = 0 (55) 
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Fic. 9 Response of airplane (B) to a 3/4 gust 
and 
yw = —O0.18 
y’ = 0.0936 


The quartic here has four real roots and the solution 
is of the form 


y = [a@sn? (wr; k) + Bb] [é sn? (wr; k) +d] (56) 
where 
k? = [(y1 — ye) (vs — ya) ]/[01 — ys) (ye vs)] (57a 
a= wie (57b 
5 = Vid (57c) 
€ = V3 — Vs (57d 
d=vyvi-—y (S7e 
and 
w = [e~°°?"/(2V 2)] V (on — ys) (v2 — yu) (S7E 


This solution (in the &, 7 plane) is plotted in Fig. 9 
and is compared with the exact solution from an isoclinic 
integration of the phase plane. The agreement in the 
overshoot is not unreasonable, and the frequency and 
damping agree well as 7 increases. 

It will be noted that this air frame possesses three 
trim points, two of which are stable (@& = 0.12, a = 
1.12); the third is unstable (@& = 0.76). The damping 
is high (¢ ~ 0.5), and the main effect of this damping 
is to make the air frame much less susceptible to any 
tendency to seek the second equilibrium point. Thus, 
in order to encounter pitch-up, much larger inputs 
are required at sea level than at high altitude. For 
example, if the same air frame were pulling up at about 


5g and encountered a —2.5g gust, it would pitch up 
violently. Here 
&% — a = —O0.18 
and 
a,’ = —1.2 


The reduced equation is 


y” + 0.488y + 1.26e7°°??y? + e738 = 0 (58) 
where 

y(0) = -1.5 

y’(0) = 1.79 








626 JOURNAL OF 





2.0 ers T ae — a. eo 
| THIRD EQUILIBRIUM POINT 
A | } 


} | € XACT(PHASE PLANE) 

















Lor 
a 
ie) T 
| 
| 
ee | = 
oO 1.0 2.0 3.0 4.0 5.0 6.0 
Tr 
Fic. 10. Response of airplane (B) toa 2(1/2)g gust while in a 
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In this case the quartic has two real roots, yi, yo, 
(y — r)? + S* [see Eq. (49)] and in terms of these 
roots the solution of the reduced equation is 
vy = jwki(v — 1) + Sy + [yeki(yr — 1) — yS] X 

cn (wr + ¢; k)} (1 —rnk+St+ 
lv. — r)ki — S]en (wr + ¢; kf (59) 
where 


kh =At VA24+1 


and the sign is chosen such that 


—(1/2)e~'*"k, < 0 
and where 
A = [s?* + (yn — 7) (ve — r)I/[sQn — ye) 
= t/t) + ky?) 


w = (ce? 2)V [S(y1 — yo) ]/(RV 1 — R?) 


g = cn , () veRi (yn = Fy te wis 4. 1.3[(y1 — 7) X 
ki + s]} 1.3[(41 — ki — s] + wok — vr) - ys} ) 


This solution has also been compared with a phase 
plane integration; while the agreement is not as good as 
the earlier example, the main features of the motion 
are recognizable and the values found by the approxi 
mation are conservative (Fig. 10) 


(5) CONCLUSION 


In the preceding pages, approximate methods were 
developed to estimate the occurrence and severity of 
of nonlinear pitch-ups in transonic swept-wing aircraft. 
The cause of such pitch-ups is usually a combination 
of tip-stall and chordwise center-of-pressure shift due 
to the appearance of shock-wave patterns on the upper 
wing surface. These phenomena, which are not dis 
cussed in detail here, lead to the appearance of one un 
stable and two stable equilibrium points on the moment 
curve. 
the first equilibrium point is applied to a trimmed air 


If a sufficiently strong initial disturbance from 


craft with such a moment curve, the resulting oscilla 
tory motion may subside about the second equilibrium 
point. This is far more likely to occur at high alti 
tudes where the damping coefficient is small and where 
the trimmed angle of attack is so large that the trimmed 
equilibrium point is close to the unstable equilibrium 


point. 
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Approximate techniques that rapidly yield approx; 
mate results and define the general nature of the ey 
pected motion have been developed to study this prob 
lem; while they do not replace analog analyzers or other 
devices for the construction of numerical solutions. 
they may serve as complements particularly useful in 
investigating the variation of flight parameters such as 
damping, undisturbed attitude, etc., analytically and 
in suggesting the most effective techniques of carrying 
out detailed computations. Numerical examples have 
been worked out to give time histories of nonlinear os 
cillations under various conditions, both at high and 
low altitudes. 


APPENDIX 
A DISCUSSION OF THE COMPUTATIONS OF THE Roots 
OF THE QUARTIC 
Consider the quartic 
x! + Bx* + Cx? + Dx + E=0 


where B, C, D, and FE are known functions of time 
If B and C are five to ten times larger than D and E, 
the quartic has two real roots, largely uninfluenced by 
D and /:, and these roots are governed by the quadratic 


x*+ Bx+C=0 


The approximation for the small roots is found by 
writing the quadratic as 


(x? ++ Bly + C’) (x? + b’n + Cc’) 0 


Expanding the product gives relations between B, B’, 
b’, ete., and letting B’ = Band C’ = C, it is found that 


c= E/C 
b’ = (CD — EB)/C? 


These expressions may be used to obtain a better ap 
proximation to B’ and C’; namely, 


B’ = B — |(CD — EB) C?| 
c’=C 
This approximation is valid whenever 
(1/C*) (E + (CD — EB) |B — |((CD — EB) /C*]|) «1 


and 
KCD — £B)/C? << D 
or if D = 0, the latter condition becomes 
(EB/C?) ~ 0 
If the roots chosen by this approximation are 1, / 
r3, and ry, then the average error in the roots is about 


Ar =~ (EF + (CD — EB) }B — |(CD EB) /C?}}) 
}CBV 1 + [2(ri2 + ro? + 753 + 7; B?); 

In the examples in the text, B = O(e*"~), C O(em" 

and = O(e“‘); hence at least two of the roots are ol 
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mation improves as the time increases. 


> 


In the second example in Section (4) for small ¢, B, 


ind C are nearly unity, D = 0, and £ is large and nega- 


tive compared to unity. This implies that it is reason 


able to assume that the roots are nearly symmetric 


with respect to the origin so that C’ = c’ = +vw7E. 
Again splitting the quartic into the quadratic and 
multiplying, 
B = b’ + B’ 
ind 
B'c! + C'b’ 0 


Thus 


This approximation is valid as long as B?/4 ~ C 
that is, for¢ < 2.0. 

Two approximations were required in the second 
example of Section (4) due to the change in the char 
acter of the motion as the trajectory is captured by the 
outer equilibrium point. 

Another approximation useful if B is of the same 
order as D and both are much larger than C and F is 
based on the existence of a real root which is approxi 


mately 
x = —B + (C/B) 
rhe small real root 1s 


x= [-—(E/D)] [1 — (CE/D?)] 


-) so that the approxi- 
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The quadratic that remains is easily found to be 


x? + [(CB (Ek, D)|x + (DB 0 

This sort of approximation to the roots has no hard 
and fast rules governing it, and, in general, once the 
approximation has been decided upon it is wise to 
multiply the factors and compare with the original 


quartic to determine its accuracy 
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The Dynamic Motion of a Missile Descending 
Through the Atmosphere 


HANS R. FRIEDRICH* ano FRANK J. DORE? 


Convair, A Division of General Dynamics Corporation 


ABSTRACT 


A method is presented for computing rapidly, yet accurately, 
the dynamic motion of a ballistic-type missile descending through 
the atmosphere. The equations of motion are separated into a 
set of ‘‘static’’ trajectory equations (zero angle of attack) and a 
set of ‘‘rotational’’ equations describing the oscillatory motion 
of the missile about its center of gravity. A transformation 
allows the rotational equations to be written in a manner anal- 
ogous to the equation for an undamped oscillating spring mass 
system with the mass equal to unity and a time variable spring 
constant. For given initial conditions this equation can be 
solved to obtain the envelope of maximum angle of attack. An 
additional transformation allows the calculation of the complete 
oscillatory motion at any time during the trajectory as a function 
of the maximum angle of attack at that time. 

This solution shows that the maximum angle of attack of a 
missile descending through the atmosphere at relatively constant 
speed is reduced even when the aerodynamic damping is neglected 


(1) NOMENCLATURE 


acceleration radial to missile flight path 


dy = 

at acceleration tangential to missile flight path 

b = reference length 

C, C2 = constants defined by Eqs. (37) and (38) 

Cp = drag coefficient, drag/qS 

Cp = induced drag coefficient due to angle of attack, 
drag/qS 

Cc, = lift coefficient, lift/gS 

Cn = pitching moment coefficient, moment /g.Sb 

Cnqg = pitch damping coefficient, dC, /d(6b/V) 

D = defined by Eq. (12) 

d = pitch damping parameter defined by Eq. (12) 

w = oscillatory frequency 

g = gravitational acceleration 

I = pitching moment of inertia about center of gravity 

Dg = Bessel function of first kind 

M = defined by Eq. (12) 

Avi = defined by Eq. (16) 

m = mass 

r = period of oscillatory motion 

q = dynamic pressure, (1/2)pV? 

R = radius of earth 

RS = reference area 

t = time 

J = velocity tangential to missile flight path 

} = Bessel function of second kind 

¥ = altitude 

PA = Bessel function defined by Eq. (33) 

a = angle between missile axis and flight path 

a = defined by Eq. (14) 

a = defined by Eq. (30) 

B = defined by Eq. (23) 

a = angle between flight path and horizontal 
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6 = lift parameter defined by Eq. (12) 
0 = angle between missile axis and horizontal 
m = moment parameter defined by Eq. (12 
é = integration variable 
p = density 
T = defined by Eq. (28) 
Subscript 


() = value associated with zero-angle-of-attack trajectory 


(2) INTRODUCTION 


— PURPOSE OF THIS PAPER is to present a simple 
yet quite accurate method of determining the 
motion of a ballistic-type missile during its descending 
flight through the atmosphere, for the case of oscillatory 
motion about the center of gravity. Knowledge of this 
motion is required for computation of the maximum 
loads, surface pressures, and temperatures, and dis 
persion of the missile from its zero-angle-of-attack flight 
path. The cause of the oscillatory motion could be 
an initial angle of attack or rotational velocity at high 
altitude or could arise from the effect of an atmospheric 
wind gust at lower adtitudes. 
cause of the disturbing condition, the designer in the 


Irrespective of the 


past has been forced to solve the equations of motion in 
a step-by-step iterative procedure due to the non 
linearity of the aerodynamic coefficients at high angles 
of attack and the time-altitude dependency of the 
dynamic pressure due to the change in both missile 
velocity and air density. 

In actual missile flights the pitching motion could 
occur in any plane, but since it is only desired to ob- 
tain the representative motion for design purposes, 
this paper assumes that pitching motion occurs in 
the plane of the missile trajectory. The additional 
assumption of a power-off, fixed controls, rotationally 
symmetric missile are not deemed restrictive for a 
ballistic-type missile. 

During the development of the German \V-2 rocket 
in World War II, H. Unger at the Technische Hoch- 
schule, Darmstadt, Germany, approached the solution 
of the problem in question in a simple way. The sim- 
plicity of Unger’s method is based on the fact that for 
many missile flight conditions the missile oscillatory 
motion about its center of gravity has negligible in- 
fluence on the trajectory followed by the missile center 
of gravity during the descending flight through the 
atmosphere. When this condition is met, the equations 
of motion describing the oscillatory motion can be 
separated from the trajectory equations. The latter 
equations can be readily solved by the usual tabular 
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methods, giving as a result the missile velocity, altitude, 
and dynamic pressure as a function of time. The 
equations describing the oscillatory motion about the 
center of gravity are then solved for the envelope of 
the maximum angle of attack. At any desired point 
along the trajectory, the angle of attack oscillation 
can be computed in detail, utilizing the pre-established 
maximum angle-of-attack envelope for initial condi 


tions 


(3) EQUATIONS OF MOTION 


The complete equations describing the motion of a 
missile descending through the atmosphere are com- 
plicated by the fact that oscillatory motion can occur 
in planes nonorthogonal to the plane of the missile re- 
entry trajectory, and with atmospheric winds affecting 
the motion. In order to reduce the problem to its 
simplest aspect, and still keep the element of practical 
application in mind, the motion is restricted to oscilla- 
tory motion in the plane of the missile trajectory and 
with no effect from atmospheric winds. 

The angular deflections, forces, and moments acting 
on the missile are shown in Fig. 1. The forces are the 
gravitational force, the aerodynamic forces due to lift 
and drag, and the inertial forces arising from the ac- 
celerations radial and tangential to the missile flight 


path, a, and a,, where 


a, = — [I cos y¥/(R+ y)] + Vy (1) 


ay = V (2 
The moments about the center of gravity are the in- 
ertial moment due to angular acceleration and the 
aerodynamic moments due to the static stability and 
the pitch damping. The set of second-order nonlinear 
differential equations describing the motion of the 


missile are then: 
—mV — CpgS — mg sin y = 0 (3) 


—-mVy + (C,/a)agS + m}[V2/(R + y)] — gf X 
cosy = 0 (4) 


629 


—I0 + Cng(Ob 1")gSb + (C, a agSb = (.)) 


6@=art ? (6 


The solution described herein depends on separating 
these equations into two sets of solutions, one describing 
the ‘‘static’’ trajectory of the center of gravity and the 
other set describing the oscillatory motion about the 
center of gravity. There are two conditions affecting 
the division of the equations. 

First, experience with the type of solutions in this 
report shows that the angle of attack generally de 
creases rapidly as the missile descends in altitude while 
the dynamic pressure increases with decreasing alti 
tude due to the rapid increase in air density. However, 
in most practical cases, the product ag has a maximum 
at some altitude before the missile strikes the ground. 
This fact entails that the induced drag due to the pitch 
angle of attack will influence the trajectory calculations 
only for an extremely short period of time and, thus, 
[This can be seen from the approxi- 
This assump 


can be neglected. 
mation that CphgS = (C,/a)a(agq)S.] 
tion is conservative for computation of the maximum 
normal acceleration because the dynamic pressure is 
larger and the duration of time in the atmosphere is 
shorter than if the actual Cp value were used. 
not amenable to 


How 
ever, there are some cases that are 
this simplification. The solution of the equations is 
then simplified by a procedure described in Section (8 
The second condition affecting the division of the 
equations concerns Eq. (4). During the time duration 
of the re-entry trajectory, the missile will generally make 
many complete oscillations. The angle of attack a, 
and with it the lift term, (C,/a)agS, will change sign 
twice during each cycle. Therefore the overall dis 
turbing influence of the oscillatory motion on the 
lateral component of the re-entry trajectory will be 
This allows the lift term to be neglected in 
Conversely, 


small. 
Eq. (4) for statie trajectory calculations. 
the term containing cos y in Eq. (4) generally has a 
negligibly small influence on the angular oscillations 
compared to the lift term. 

The system of Eqs. (3) through (6) therefore reduces to 
mg sin yy = O (7) 


—m |" Co,GuS 


—mVoyo + mj (Vo? (R + y)] — gj cos yo= 0 (S 


which describe the static trajectory of the missile, and 


= Voy + (Cc, a) ag S=0 (Q) 
—16 + Crg(6b/ Vo)qoSb + (Cy, a)agSh = 0 (10) 
6@=act 7 (11) 


for the angular motion of the missile about its center 
of gravity. 

A tabular solution of the static trajectory equations 
is relatively simple and straightforward. For the re- 
mainder of this paper it will be assumed that the static 
trajectory solution is available, hence that the time 
dependency of go, Vo, m, and J is known. 
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(4) SOLUTION FOR THE ENVELOPE OF MAXIMUM 


ANGLE OF ATTACK 


The following definitions* are used throughout the 


remainder of the paper. 


56 = (C,/a) [(qoS)/(mVo) | 

d Cing{(GoSb7)/ (1 Vo) | 

oo = (Cm/ a) (GoSb/T) (12) 
D=d+6 


M=ypt+d+6 


Using these definitions allows the combination of Eqs. 
(9), (10), and (11) to give 


a+ Da+ Ma = 0 (13) 


The following deduction 1s based on a method described 
Introducing the transformation 


Sl Pe 1g 
= (14) 


in reference |. 


a(t) = a(tje 
reduces Eq. (13) to 
a+ Ma =0 (15) 
where 


M = M —- (1/2)D — (1/4)D? (16) 


Eq. (15) corresponds to the equation for an undamped 
oscillating spring mass system with the mass equal to 
unity and a time variable spring constant, .1/(f). 
Multiplying Eq. (15) by @ and adding J7a@*/2 to 
both sides of the equation allows it to be written as 


(d/dt) [(@?/2) + (Ma/2)|] = Ma?/2 (17) 


This equation says that the time rate of change of the 
total energy of the “‘spring mass’’ system- namely, 
the sum of the kinetic and potential energy- is equal 
to fa?/2. 

For the maximum and minimum points, @nax, of the 
a curve the kinetic energy is zero because @ = 0. The 
total energy then consists of potential energy only 
namely, MJ @max*/2. By replacing the total energy 
(a@?/2 + Ma*/2) on the left side of Eq. (17) by a fic- 
titious potential energy J7@max”/ 2 for the same ‘‘spring 
constant”’ 17, the differential equation of another curve 
is found which can be considered as the locus for the 
maximum and minimum points of the &@ curve. This 
fictitious potential energy represents the true poten- 
tial energy of the @ curve at the maximum and mini- 
mum points. Eq. (17) can then be represented by 


(d/dt) (Mamax?/2) = Ma?/2 (18) 

or 
@nax/ @max = —(1/2) (A/M) [1 — (&/amax)?] (19) 
The curve described by this differential equation still 


has a big disadvantage. It is not a smooth connection 


* The time dependency of the functions is noted in the re- 
mainder of this report only when required for clarity 
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of the maximum and minimum points of the & cury; 
Because of the term @/@n.x” in Eq. (19), the locus 
CUrVE Ayax 1S Slightly oscillatory, displaying double thy 
frequency of the &curve. 

For ballistic-type missiles, the assumption generally 
can be made that the functions d, 6, w and thus 17 do not 
show extremely rapid changes over the time interya] 
required to complete half a cycle of oscillatory motion 
Therefore, by only slightly affecting the location of th: 
maximum and minimum points of the & curve it is pos 
sible to assume the & curve approximated from half 
wave to half wave by a sine oscillation as 


Q@/Qmax = Sin wi 20) 


From Eqs. (19) and (20) then follows 
max/@max = —(1/4) (7/17) (1 + cos 2ef) (21 
This last equation clearly shows the wave form oi 
the locus curve and that it makes one full cycle for 
each half cycle of the @ curve. It follows that the 
mean value for each half cycle of the & curve is then 


l 1 /I 4 
- “ | cos 2wtd(wt) = sin 2w/ 0) 
T/oJSO0 T 2 2 0 


The wave curve @,,, can then be smoothed by the 
following excellent approximation 


Gmax/ Bmax = —(1/4) (M/M) (2: 


From this equation follows the equation for the 
envelope, at any time ¢ > 0, of the maximum and mini 
mum points of the & oscillations described by Eq. (15 


namely, 
Amax(t) = Amax(O) ¥ [I7(0)]/ (AMD | (24 


In Eq. (24) the time ¢ = 0 is defined such that &,.x(0 
coincides with the first maximum or minimum of the 
& curve after it has become possible to approximate the 
& curve from half wave to half wave by the sinusoidal 
motion indicated in Eq. (20). This condition is illus 
trated in Fig. (2a). 

The equation for the envelope of the maximum and 
minimum points of the angle of attack oscillations 1s 
found by combining Eqs. (14) and (24) to get 


l 
! pitydt 
So D(é) dé 


4 _ = 9 
Qmax(t) = @max(O) V (M(0)]/[M(b]e ° 


(25 


but from Eq. (14) and the above definition of the time 
t = 0 follows @nax(0) = [a(O)] max = a(0). Therefore, 
the general equation for ay, (¢) is 


Si D(é)dé 
2 JV (96 


Qmax(t) = a(O) V {AN(0)]/ [Je 2 


(5) MoTrion DuRING A CYCLE 


After the envelopes of the angle of attack oscillations 
have been found according to the method in the pre 
ceding section, it is possible to calculate the representa- 
tive angle of attack oscillations, piecewise, in a closed 
form at various points of interest. 
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The solution for the maximum angle of attack enve- 
lope of the preceding section is general in that nonlinear 
aerodynamic coefficients can be handled. In this sec 
tion the solution for the oscillatory motion during a cycle 
is more restricted in that it is assumed that the param 
eter .J/(t) can be approximated with sufficient accuracy 
over a quarter cycle by the relation 


M(t) = M(t) (27 


where 
r = ef" (28) 
and where f) is the time at the start of the quarter 
eyele. [Eq. (27) is exact for a missile with linear 
aerodynamic coefficients descending at constant ve- 
locity in the isothermal atmosphere.] 98 is obtained 
from Eqs. (27) and (28) as 
8 = [1/(t — t)] In, M(t) / 1 (to) (29) 
For the conditions where Eqs. (27) and (28) are appli- 
cable, the parameter is relatively independent of time. 
Introducing the transformation 


a(t) =al(r (30) 
allows Eq. (15) to be written as 


r|d2a(r)/dr?] + r[da(r) /dr] + 
[\7(t)) B2]a(r) = 0 (31) 


which has as the solution” 


a(r) = Zo|(2/8) V M(to)r] (32 
where 
Lol CiJo( ) + Co ¥ol (33 


J) is the Bessel function of the first kind of order zero, 
and J), is the Bessel function of the second kind of order 
zero, both of argument |(2 8) V 17 (t))7]. 

Eqs. (14), (27), (28), and (30) are used to transform 
Eq. (32) into an expression for a giving 


45 Deeg - 


a(t) = Zo[(2/B) V A(t) Je (34) 
The values of & and & are obtained from Eq. (34) as 
a = }—(1/2)DZ,[(2’/8) VM] —- 
Sa 
VM Z,{(2/8) V AT] e afo Pwat (35) 
& = }[(1/4)D? — (1/2) D — MJZ,[(2/8) V M1] + 


a 5 fF Diedé 
(2/8) V M)j e 2 So (36) 


DVM Z, 





where Z; is the Bessel function (or cylinder function) of 
order one. The values C; and C, of Eq. (33) are deter- 
mined by assuming a time f) such that a(f)) = O, there 
fore a(to) = Qmax(to), Where @max(to) is taken from the 
maximum angle of attack envelope. Substituting 
these conditions into Eqs. (33), (34), and (35) and solv 


ing for C, and C. gives 


C, = Amax(to) — V Mite) Vi[(2 8B) V M(t) ] + C1. 2)D (to) Yo[(28) V M1 (te) | / 
Co = = = = = +4 
V M(to) | J[(2'8) V M(t) ] ¥i((28) V W(te)] — Aif(28) VY AM (te)] Yo((2/8) V AM (te) J 
C Amax(to) — | V M(t) Af(2/8) V M(t) |] + 1 2)D (te) Jo (28) V A (to) | / 
2 = iz x = = 3S 
V M(t) | —J[(2'8) V M(to)] ¥il(28) V M(t) ] & AL(2 8B) VY AM (to) |Vo[(2 8) VAT (ty) JN 
It should be noted that a@4x(¢)) appears in both C; and t = ty; cos [(2/8) V W(t)] = (39 
(:, therefore in Z) and Z;.. From the form of Eqs. (34), 
39), and (36), this means that a(t), a(t), and &(t) are at 
6 seemnghe ase ei t= t + P/4; cos [(2'8)V. Mm + P41] =0 (40 


Once the values of a, a, and & are known, the values 
of 6, 6, and @ are readily obtained by the application of 


Eqs. 6) and (12). 


(6) FREQUENCY OF THE MOTION 


In the preceding section it was shown that the motion 
during a cycle could be approximated by Bessel func 
tions with the argument [(2 8)V/J/(t)]. For large 
values of the argument, the Bessel function can be 
written in terms of sine and cosine functions with the 
same argument. This means that the oscillatory por- 
tion of the motion can be represented by sinusoidal 
motion. For example, if the angle of attack is a maxi 
mum at ¢ = fe, a quarter of a period later it will be zero, 
or writing this in terms of the cosine of the argument, 
at 


But since a quarter period corresponds to 7 2 radians, 
the argument of Eq. (40) can be written as 


(2/B)V M(t) + (2/2) = (2/8)V M(t + P/4 4] 


But from the approximation for 17(¢) given by Eqs. 
(27) and (28), the period can be obtained from Eq. (41) 


as 
P $8) In, } {VW W(t) + 28/4]? [10(%)]! (42) 


with 8 given by Eq. (29 

Since the period is time dependent it is desirable to 
be able to plot the envelope of the period as a function 
of time, ¢. The derivation for the period given above is 
still correct if the parameter @ is evaluated as 
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M(t) : Sa an 
8 = lm ni. = (43) 
tt L(t — bo) M (to) 
emax ENVELOPE 
With this definition of 8 it is then correct to speak of 4 


P(t) as the instantaneous value of the period at time f, 
analogous tO Qmax(t) being the instantaneous value of 




















Q@max at time ¢. The oscillatory frequency is given by 2 
a “ - a = 
the usual definition \ ~aiee —— 
—" 
+ an > IM 0 
w = 2n/f (+4) —* _—— +, TIME INCREASING 
a —= ALTITUDE DECREASING 
— 
(7) APPROXIMATE SOLUTIONS FOR THE DyNAMIC REGION OF APPROXIMATE rh 
Mowrson SINUSOIDAL MOTION ie 
éyax ENVELOPE — > torqt 
The Bessel solutions for the dynamic motion a = a \ b tions 
; a ae ee TV i ; 
a(t), a(t), and a&(t)—are cumbersome to evaluate. +i} — ously 
rm1* “~. . . - \ - t ing ¢ 
rhis difficulty is overcome by taking advantage of two _ ~ a 6 
ie ae : ig ss =e Si famil 
conditions indicated by the Bessel solution: first, gs wit i 
a(t), a(t), and &@(¢) are proportional to anax(t), and ad 
second, the oscillations can be approximated by sinu- ZY = tory 
soidal motion with a time dependent frequency, w(t). diyax ENVELOPE P \ 
These statements lead to the type of equation ; “ a c 
- ao Sees \ 
MAX eaten 7. j 
- ANZ a es 
a(t) = Qmax(t) Cos [(w(t)t] (45) - 
ME +=0 a 
~ / 
Differentiating this equation for a(t) and a(t) gives ™ ’ 

a VR 
~ Vos , i ifs, ; = ; ; —_ re 
a(t) = — max, Sin [w(t)t]j [w(t) aie w(t)t] 45 Fic. 2. The dynamic motion of a missile descending through weal 

Gmax) COS[w(t)t]} (46) the atmosphere iden 
} co p 
&(t) = —aimaxycos [w(t)t]} [w(t) + a(At]? — the aerodynamic damping term is small, then certain | — abot 
Qmax}sin [w(t)t]} [Qe(t) + w(t)t] — useful approximations can be made to simplify the whe 
ae ; calculations considerably. Generally speaking, it may | i 
Pémax) sin [w(t)t]}} [w(t) + w(e)t] + . - : y ing, may | axis 
te _ be found that | func 
Q@max COS [w(t]; (47) 
: " . ~ | sup 
u>b6+ ds — (1/2)D — (1/4)D 51) PI 
where RG 
so that 
Qmax = [damax(t) | dt / It 
‘ (48) ao a — the 
o 9 9 j = (02) | 
Qmax > [d?Qmax(t) | dt?\ B the 
pa . on : is an extremely good approximation. Furthermore, gene 
The maximum values of a(t) and &(¢) are usually easier . : ‘ 
: ae : the aerodynamic damping builds up slowly so that the 
and more useful to work with. These are 3: : 
5 Se Dé) dé w it 
; , 290 = — ie 
[ a(t) Jmax = Qmax(t)w(t) (49) é <4 Dd) sign 
and for a sizable portion of the re-entry trajectory. 
° , or . With these assumptions Eq. (26) can be approximated 
[a&(t) max = Omax(t) [w(t)]? + [d2amax(t)]/dt?} (50) : I | PI 1 
by 
It is seen that the equations for [@(t) max and [@(t]max i type 
“ete : e ; : sea \ e : 
are similar to the values obtained for simple harmonic Amax(t) ~ a(O) ¥ [H(0)]/ [a] of wit! 
motion except for the last term of Eq. (50). In prac- “—™.: er . : 3 ; 5 ticle 
; : Pp - /. \ I rhis equation is readily solved by iteration. Eq. (54 Lay 
tice this particular term is usually of small magnitude. — gi . - vers 
: is interesting because it clearly shows that the maxi- wm 
» ss . OV 
: : mum angle of attack is usually reduced even when the 
(S) APPROXIMATE SOLUTIONS NEGLECTING ; ; , age tanc 
aerodynamic damping is neglected. This occurs be 
AERODYNAMIC DAMPING é ;: , . 
cause the increase in dynamic pressure due to the 
One of the original suppositions of this paper was density change with decreased altitude effectively 
that the solutions are particularly applicable to high- amounts to an increased ‘‘spring constant”’ in the mo- whe 
speed missiles in their descending flight through the ment equation, requiring a decrease in angle of attack righ 
atmosphere. Now for some types of missiles it is pos- (i.e., ‘spring deflection’) to maintain the same kinetic g 1S 
sible that the ratio of the aerodynamic damping term energy of the system. Obviously the opposite situation the 
to the static stability term of Eq. (5) is small, in fact, can occur (pitch oscillation divergence) if the reduction 
much less than for airplanes or missiles designed to in velocity due to a large amount of drag is more than 7 


maneuver at transonic or low supersonic speeds. If (Continued on page 638) 
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Some General Comparisons Between the 
Vibratory and Conventional Rate Gyro 


JOHN B. CHATTERTON* 


Sperry Gyroscope Company Division of Sperry Rand Corporation 


ABSTRACT 


The new with 


conventional elastically restrained 


vibratory rate gyro is compared analytically 


rate gyro. Developed 


transfer func 


thie 


torques, sensitivity bandwidth factors, output 


tions, and some of the inherent errors are developed simultane- 
ously in a simple unified manner to give a theoretical understand 
ing of the vibratory gyro and its characteristics relative to the 
familiar conventional rate gyro. The discussion is of a general 
nature and does not attempt to compare performance data of 


practical instruments. The conclusion is reached that the vibra- 


tory rate gyro has many relative advantages 


INTRODUCTION 


R* ENT DEVELOPMENTS in the field of gyroscopics 
brought forth the vibratory rate gyro 
, which has the obvious advantage of no moving 


have 
VRG 
wearable parts. 
identical masses, electrically maintained in opposed 
co-planner vibration, produce an alternating torque 


Its usual configuration is such that two 


about the symmetrical axis in the plane of vibration 
when angular inertial velocity is applied about the 
axis The device is comparable in its measurement 
function to the conventional elastically restrained (and 
supported) single degree of freedom gyro or rate gyro 
RG 

It is the purpose of this paper to compare some of 
the more interesting and important characteristics of 
the two devices. 
general discussion of inherent properties rather than 


The comparison is held to a very 


the examination of characteristics of working models. 
It is believed that what follows is considerably more 


significant. 


DEVELOPED TORQUES 


In the derivation of the torques developed for each 
type of gyro, the similarities are best shown by starting 
with the familiar vector equation of motion for a par- 
ticle in spherical coordinates. The radial and trans- 
verse components of velocity of a particle of mass, m, 
moving in a secondary inertial system, x, y, 2, at a dis- 
tance r from the origin is (see Fig. 1): 


v= tr + r69 + rd sin 06 (1) 
where * is the unit vector along 7, @ isa unit vector at 
right angles to 7 and in the direction of increasing 6, and 
¢ 1s a unit vector at right angles to the *@ plane and in 
the direction of increasing ¢. The time derivative of 
Received June 3, 1954. 

* Senior Engineer, Flight Instrument Engineering Department 
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this velocity’ gives the vectorial acceleration of the 


particle relative to the inertial system 


a= ly — r(0)- — rio sin? @|7 4 
[276 + 76 - r\o sin 8 cos 6\6 — 
[2% sin 06 + 2r@¢ cos 6 + r@ sin A\o (2 


The typical rate gyro configuration consists of a 
rigid rotor mounted in an instrument frame and oriented 
with respect to a spin, input, and output axis as shown 
in Fig. 2. With the instrument frame axes positioned 
with respect to the inertial system as shown, and rotat- 
ing in it with velocity 2, Eq. (2) is valid for each mass 
particle of the rotor. The reaction torque output of 
any particle will be due only to the @ component of ac 
celeration, therefore, with 7 of constant length, this 
incremental torque amplitude becomes, with the direc 


tion of positive output, 


t = m[2r6d cos 0 + r@ sin 6|r cos 8 2 
With @ = © the inertial angular velocity input, with 
6 = a,» the rotor velocity, and considering @ 0 for 


the moment, Eq. (3) becomes 
t 2mr-wr2 cos~ 6 (4 
If the total rotor mass, .1/, is considered to be distributed 
along the circle at r, 
m (.\/ /22)dé (5) 


and the total output torque becomes 


r ‘wre MQ 


eon 
} cos? 6 da (6) 
1 aa 


= rwrMQ (7 


Tre 


or 
Tre 


which, of course, is simply /wr2, where / is the inertia 
of the rotor about the spin axis. 

The vibratory rate gyro configuration is somewhat 
similar and is shown in Fig. 3 with rotation in the in- 
ertial frame as before. The instrument frame coor- 
dinate axes can be designated here as the input, vibra- 
tion, and transverse axis. The mass vibrates about the 
vibration axis in the input-transverse plane and sym- 
The output 
To avoid a 


metrically with respect to the input axis. 
axis could be taken along the input axis. 
double frequency output torque, however, the output 
axis is taken parallel to the input at a transverse dis- 
tance R. The output torque will again be due to only 
the @ component of the acceleration of Eq. (2). The 
magnitude of the incremental reaction torque output 
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in the direction of negative input becomes, with r con- 


stant in length, 
t = m[2r6¢ cos 6 + r@ sin 6| (R — r sin @) (S) 
If the motion of 77 1s such that 
= AO sin wrl (9) 


representing a mass vibration of peak amplitude A@ and 
angular frequency wy, and with 6 = 2 and ¢@ = O as 


before, Eq. (S) becomes 


f = m[2rQapA0 cos wrt cos (AP sin wrt) | X 
[R — r sin (AO sin wet)| (10) 


The total vibration mass / can be lumped at m and 
assuming A@ small, the total output torque can be 
written approximately, 


7\ RG ye Rr Abc, 2 cos wpl (7d) 


Since the input angular acceleration @ was considered 
to be zero, Eq. (11) will still be valid for 2 applied 
about the z’ or output axis. 


OutTPeut MOTIONS 


It is the usual practice to have these torques produce 
a torsional motion of a restraining system with respect 
to the instrument frame and, with some appropriate 
transducer, produce an output voltage or signal pro- 
portional to the motion. Typical configurations are 
shown in Figs. 4 and 5. 

The RG must utilize a differential displacement sen- 
sitive transducer in order to detect differential applied 
rate; hence, the displacement of a spring mass system 
driven by the torque of Eq. (7) will be calculated. 
With /) the total moment of inertia about the output 
axis, k, the torsional spring constant of the torsion 
system, and C the mechanical damping resistance, the 
displacement, Dy, for the steady state case when 2 = 


1” cos wl, becomes 


Dac = Tre ljeljols + C + (k/jw)}} (12) 
With the natural angular frequency, 
w = VR I, (13) 
and the Q of the system 
QO = (wlo/C) = 1/[2(C/C,.)] (14) 


where C, is the critical damping resistance, Eq. (12) 
becomes 
Dre 
Derg = cs . (15) 
R{1 — (w wy) ~ fe Kl Q) (w wo) | 


Substituting the torque of Eq. (7) and converting the 
imaginary expression 


Dre = (PwrMV kg) cos (wt — 8) (16) 
where 
8 = tan! [way/Q(wy? — w?) (17) 


and 


g = Vl — (w/w)?}? + [1 /Q) (was) | IS 


For a practical design the Q is usually chosen near ().(2 
as a compromise for best response time and _ transient 
overshoot. For this value the approximate rate trans 


fer function becomes, from Eq. (15), 


Gro(jo) = Tre kl + jw/eo) | 19 


resulting in a maximum phase shift of 1SO° lagging as 
would be expected in a second-order system. 

The VRG, because of the higher frequency carrier 
centering of torque, can utilize a velocity sensitive 
nondifferential transducer and, hence, the velocity of 
spring mass system driven by the torque of Eq 


will be shown. 


With 2 = I cos wl as before, 
wre = Tra |jolo + C + (Rk jw) | 20 
Using Eqs. (13) and (14) with wy = ow») and assuming 
g Eq 


arithmetic symmetry of the response near resonance so 
that 


| co Oy TT W 2(w ‘wy A 21 


| (ax T Ww wy | 


the solution for Eq. (20), assuming w < w), can be shown 


to be 


MRr Abo V 
CV1+ OK? 


UVAG 


leos|(ws + w)f — a] + cos |(w& w)i +a 

where, a = tan 'QA. Eq. (22) shows the double side 

band characteristic of the suppressed-carrier motion, 
and rewriting, 

2A Rr Ade V 

cos (wl — a@) COS wl a 


CV1+ OK? 


UVRG 


The rate transfer function of interest would be that of 
the velocity output after a perfect demodulation or 
phase detection operation has been accomplished. 
This could be written immediately from Eq. (25) after 
the cos wf factor is dropped and w is assumed small to 
be consistent with the condition on Eq. (22). This 
condition could then be changed to give an approxt- 
mate output velocity for large w, and complete the 
transfer function approximation. However, a simple 
and readily visualized method, similar to a.c. servo 
techniques, can be used to obtain the entire output re- 
sponse. If the torsion system amplitude response char- 
acteristic is assumed to be a perfect even function about 
its natural frequency, w, and the phase, a, a perfect odd 
function about w,, approaching 90° lead to the left of 
w, for lower frequencies and 90° lag to the right, then 
for an arbitrary suppressed-carrier input, centered at 


Wiis 


Input = A cos wf cos wolf 
(A/2) [cos(w + w)t + cos (wy — wt] (24 


and, with points AA on the amplitude curve corre- 
sponding to frequencies w) + w and wo — w, the output 


is 
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Output = (A/2)AA }cos [wot + (wt — a)| + 
cos [wot — (wt — a)}} 


= AAA cos (wl — a@) COS al (25) 


Hence, the demodulated output is identical to the right 
or higher frequency symmetrical half of the torsion 
response curve with the value of w = 0 starting now at 
a. This is clarified in Fig. 6. When these character- 
istics are redrawn on semilog coordinates, the approxi- 
mate attenuation asymptotes can be drawn immedi- 
ately. This is shown in Fig. 7 along with the phase 
characteristics for each asymptote. The solid lines 
show the true value and the dashed the approximate. 
Thus, from Eq. (23), and with the restriction on w as 
seen from Fig. 7, the transfer function can be written 
for the demodulated output 


: : Tyre/COS aot : wy : 
Grre(jwo) = , for w< (26) 
1 + 20](w/wo) 20 
and 
; : 1 yRrG/ COS wot - ww oa 
Grre( jw) = = for w > (27) 
207] (w/o) 20 


Thus, the maximum phase shift possible in the I’RG 
is 90° lagging and can be approximated by a first-order 
rather than second-order system. 


SENSITIVITY BANDWIDTH FACTORS 


In any transducer it is generally desirable to have as 
large a sensitivity (using sensitivity as a gain factor) 
as possible. Furthermore, if the transducer or instru- 
ment is required to transmit information, it is desirable 
to have as wide a frequency response as_ possible. 
Therefore, the possibility exists that the product of 
sensitivity and bandwidth for each gyro could be a use- 
ful figure of merit. 

The significant bandwidth for the RG is simply 


BWreg = o/22 (28) 


and the sensitivity, from Eq. (16), per unit steady rate 
input, 


Sre = (rwrM)/k (29) 
or, in terms of the spin and output axis inertias, 
Srq = (wR/a*) (I/Lo) (30) 
The sensitivity-bandwidth product is then 
SBWre = (wr/27) (1/1) (31) 
For the |’RG, from Eq. (23), the sensitivity is 
Sveg = (2MRrAdwp)/C (32) 


which is a function of the damping unlike the RG. 
Setting AR = rAd, the mass transverse displacement, 
and, comparable to the spin inertia of the RG, a vibra- 
tory axis inertia of 7 = .VR?, the sensitivity becomes, 
using Eq. (14), 
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Si = (2ARQ R) (J is) 33 


The bandwidth, considering the double sideband nature 
of force, as was seen in Fig. 7, 


BWrre = ‘4drQ 34 
and the sensitivity-bandwidth product becomes 
SBW BG = (ARw; 2rR) (1/1) 35 


The result for the RG, unlike that for the RG, js 
seen to be constant for any particular design. The 
factor for the RG depends on the inverse bandwidth 
(a high price in sensitivity is paid for an increased band- 
width) and is therefore not a particularly useful figure 
of merit except for this comparison. As an example of 
the numerical values of these results, take the I’RG 
with a vibration amplitude of 0.002 in. at 2,000 cycles 
per sec., an R of 0.25 in. and an inertia ratio of 0.3; for 
the RG, take a rotor speed of 12,000 r.p.m. and an 
inertia ratio of unity. The two SAW factors will be 
equal for an RG bandwidth of 6.6 cycles per sec., the 
SBW factor for the VRG being higher for RG band- 
widths larger than this value. For equal pickoff 
transducer gains (volts per unit angular displacement 
and volts per unit angular velocity), these figures 
would represent the relative electrical outputs for the 
same applied angular velocity inputs. 


SUBMINIATURIZATION 


The sensitivity-bandwidth factors can be used as 
general design equations for possible subminiaturization. 
As the configuration for either device is scaled down in 
physical size, the inertia ratios will remain constant. 
With constant rotor speed for the RG, the sensitivity 
and SBW factor will be independent of physical di- 
mensions if the bandwidth is held constant. The 
bandwidth is likely to increase with subminiaturiza- 
tion, however, and, depending on the application, might 
result in too low a sensitivity. The VRG, on the 
other hand, will have the product ARw, nearly inde- 
pendent of physical size and a decrease in R will in- 
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crease the SBIV factor. Scaling down the output 
pickoff will lower the overall sensitivity in each case. 
Subminiaturization of the pickoff will not affect the 
performance of the VRG through a resolution loss, 
however, since it utilizes a nondifferential pickoff. 


[his will not be the case for the RG. 


ERRORS 


It is not necessary to point out that the output 
motions discussed are for the ideal case, where output 
js linear with applied inertial rate. There is probably 
an uncountable number of deviations from the ideal, 
quantities which can be added or superposed on the 
ideal and quantities which become factors of the ideal. 
These deviations can perhaps be divided into two 
groups uncertainties, those deviations that are diffi- 
cult to discuss analytically in the general case; and er- 
rors, those that more or less may be described quan- 
titatively. In general, the uncertainties limit resolu- 
tion. In comparing two different types of instruments 
itis therefore mandatory that a discussion of the largest 
or most important of these deviations be included. 

In developing Eqs. (7) and (11), the applied angular 
acceleration about the input axis, @, was assumed to be 
zero. Considering this component now, for the case 
of the RG from Eq. (5) 


f, = mr*g sin 6 cos 6 36) 


and the total component becomes 


> 


Mro (2 
: | sin 6 cos 6 dé = O 


TRG = 2 


=i 


For accelerations about the output axis, however, there 


is obviously 
an en 
Trea = Lodout (38) 


which would be superposed on the developed torque 
due to an input velocity and represents an acceleration 
error. 

For the VRG, the acceleration component becomes, 


from Eq. (S 
ts = mr@ sin 0(R — r sin 8) (39) 


and using Eq. (9) with A@ small, the total torque be- 


comes approximately 
Tyres = MRrAd¢ sin wrt (40) 


In addition to this component, transfering the input 
axis to the s’ or output axis will cause an angular ac- 
celeration component of torque, similar to Eq. (38), 


Tyres = lod (41) 


The two torques of Eqs. (40) and (41) represent an ac- 
celeration error for the VRG about its input axis. It 
is seen, however, that the torque of Eq. (40) is out of 
time phase with the primary torque of Eq. (11). It is 
therefore possible to eliminate this quadrature term 
When the output signal is subjected to a demodulation 
process. The component of Eq. (41) is not of the car- 


rier frequency and would ordinarily be attenuated by 
the torsion system response characteristic. 

More generally, the acceleration errors discussed 
above become transient errors when the applied input 
is not considered to be a steady sinusoidal. For a step 
function of velocity applied simultaneously about the 
input and output axes of the RG, the output response 
will be the response of the output torsion system for 
torques of step and impulse functions respectively. 
In fact, the latter output amplitude can be made rela- 
tively low by proper design. For the step function of 
velocity applied to the input axis of the VRG, the fore- 
ing function of the torsion system will be a step function 
of carrier torque and a unit impulse due to Eq. (41). 
Since the Q of the torsion system is high, an intolerable 
overshoot would result. It is therefore mandatory that 
either the Q be lowered or that the instrument frame be 
isolated from the vehicle and offer low transmissibility 
to frequencies near wr. Since wy is usually large com- 
pared to the BW, the latter solution is very practical. 

Another source of possible error in each device 1s 
the effect of output motion on the applied input. 
For the case of the VRG, the output velocity of Eq. 
(23), vvre, is applied to the input like feedback. The 
corrected output velocity is then 


UVRG UvRG 
ila = tyre | 1 — 
1+ y RG Q Q) 


which reduces approximately to vyrq since the second 
term in the denominator (although indeterminate at 
times) is always less than unity. The displacement 
output Dre of the RG, Eq. (16), causes a displacement 
between the effective input axis and the prescribed in- 
put axis of the frame. The corrected output displace- 


ment becomes 


dD, = Dre cos Dre = Dre {1 —_ (Dre 2I)+.... 


which reduces to Deg for small displacements. Com- 
paring the expansions of Eq. (42) and (43), both output 
motions appear to be nonlinear. The VRG output, 
however, is seen to be linear if only the carrier fre- 
quency term is used. This could be practically ac- 
complished with the use of a full-wave demodulation 
process on the output. The output motion of the 
RG will have another effect in that applied rotation 
about the frame spin axis will cause an additional out- 
put displacement proportional to the sine of the total 
displacement and will result in a cross-coupling error. 
Another error of the RG, which has nothing com- 
parable contained in the VRG, is the error due to me- 
chanical hysteresis of the output torsion system. 
Essentially the same system is used for the VRG 
but because of the carrier-type of torque, hysteresis 
produces only a negligibly small phase distortion. 
Because of the errors of cross-coupling and hystere- 
sis, an upper limit on allowable output displacements 
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for the RG is required. Consequently, the VRG, with 
no upper limit on output motion other than the fatigue 
endurance limit, has a practical linear amplitude range 
of perhaps 1,000 times that of the RG. 


UNCERTAINTIES 


Since the inertial rate that is desired to be measured 
is a differential quantity (clockwise and counterclock- 
wise directions must be distinguished) a device to meas- 
ure this must be of the differential type. Any prac- 
tical differential measuring device will have a nonzero or 
minimum null with zero applied input. Since the cause 
of these nulls could never remain entirely constant, 
in order that compensation or bucking could be totally 
effective, a source of a null will also be a source of vari- 
ation of the null or an uncertainty. 

The RG has essentially three sources of null and un- 
certainty. These are due to the static unbalance of the 
output system about the output axis when acted upon 
by linear accelerations such as gravity, the noise pro- 
duced by the dynamic unbalance of the rotor, and the 
displacement pickoff which must be of the differential 
type. 

The VRG has a nondifferential pickoff and a mono- 
which produces no_ noise. 


chromatic motion 


There are nulls, however, due to the dynamic unbalance 


mass 


of the double-ended vibrating system and to a possible 
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static unbalance when acted upon by linear acceler 
ations. It has been found experimentally, however 
that this later null is apparently a function of the first 
and follows the dynamic unbalance as a second-order 
effect. 
essentially one source of null and uncertainty. 


Thus it might be said that the VRG contains 


CONCLUSIONS 


There is no doubt that the VRG has inherent features 
relative to the RG which make it outstanding in many 
areas. In fact, an audit of absolute advantages and 
disadvantages for each device would certainly favor 
the VRG. 
usefulness for application, especially today’s applica- 


If the audit were weighted according to 
tions, the superiority is not as great. The fact, how- 
ever, that the VRG is a new device with a lower number 
of causes of uncertainty should make it relatively more 
susceptible to future development. 
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sufficient to offset the increased density at lower alti- 
tudes. 

Approximate expressions can also be obtained for 
fae(t) |inax, [&(t) max, and w(t) with the same restrictions 


as for Eq. (54). The resulting equations are 


(a(t) Jax SE [&(O) Jnaxy (u(t) ]/ [u(0) |} ws (55) 
[&(t) Jax SE [&(O) Jrnax} [e(t) ]/ [u(O) ]} 7 (56) 
w(t) & w(O) } [u(t)]/[u(O) ]} (57) 


The use of Eq. (54) allows a rapid calculation of the 
maximum possible @nax envelope, which can then be 
used to establish the magnitude of D(t) and J/(t) re- 
quired for the exact calculation of ama, with the con- 
sideration of aerodynamic damping. In addition, the 
effect on the ‘‘static’’ trajectory of the increased drag 
coefficient due to angle of attack effects can also be 
readily computed. 


(9) CONCLUSIONS 


This paper presents a method for computing the 
dynamic motion of a ballistic-type missile descending 


through the atmosphere. The equations of motion 
are separated into a set of ‘‘static’’ trajectory equations 
(zero angle of attack) and a set of ‘“‘rotational’’ equa- 
tions describing the oscillatory motion of the missile 
about its center of gravity. These equations are 
solved separately to obtain the 


and the envelope of the maximum angle of attack. At 


‘static’ trajectory 
any desired point along the trajectory, the angle of 
attack oscillations can be computed in detail, utilizing 
the pre-established maximum angle-of-attack envelope 
for initial conditions. The type of information re- 
sulting from this calculating procedure is shown in 


Fig. 2. 
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Criteria for Longitudinal Stability 
Requirements 


JAMES L. DECKER®* 
The Glenn 1. Martin Company 


ABSTRACT 


This paper proposes a new set of requirements for the longi 
tudinal stability characteristics of aircraft. The factors in- 
fuencing the flying qualities of piloted aircraft have been an 
lyzed, and the proposed requirements are based primarily on 
stick-free dynamics and control force characteristics It is 
shown that the stick-fixed neutral point is not always a valid 
riterion for longitudinal stability 

The proposed requirements are listed below 

The aft center-of-gravity limit (1) shall not be behind the 
g. location at which the stick force-per-g is less than a specified 
value, depending upon airplane type and limit-load factor, and 
2) shall be forward of the stick-fixed short-period divergence 
point 

The short-period dynamic stability characteristics shall be 
such that: 1) the stick-free short-period airplane oscillation 
lamps to one-half amplitude within one cycle and (2) the stick 


free short-period airplane oscillation damps to one-half ampli 


tude within 1! sec 


The latter requirements may also limit the aft center of gravity 


loading 


rhe long-period stick-free dynamics may be unstable—except 


in cruising and climbing flight—but the motion shall not double 


In cruising and climbing flight, 
ositive-force stability must exist so that push and pull forces 
re required to stabilize at speeds above and below trim, respec 


tively. This force gradient should be of sufficient magnitude to 
maintain essentially the trim speed after a sharp control dis 


placement and release 


SYMBOLS AND DEFINITIONS 


D = drag, lb 

‘ = lift, lb 

Vi = pitching moment, ft.lb 

F = control force, Ib 

‘i = period of oscillatory motion, sec 

.) = wing area, sq. ft 

1 = thrust, ib 

V = velocity, ft. per sec 

W = airplane gross weight, Ib 

( = mean aerodynamic chord of wing, ft. (MAC 


= gravity acceleration, ft. per sec.” 


n = load factor in g's 
= dynamic pressure, lbs. per sq.ft 
I, = pitching moment of inertia, slug-ft.? 
a = angle of attack, radians 
6 = pitch angle, radians 
6, = elevator deflection, radians 
C, = lift coefficient, L/Sg 
Cy = drag coefficient, D/Sg 
Cn = pitching-moment coefficient, .11/S.q 
Cua = the slope of the elevator-hinge moment-coefficient 


curve with tail angle of attack at constant elevator 


deflection 
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Cys = the slope of the elevator-hinge moment-coefficient 
curve with elevator deflection at a constant tail 


angle of attack 


Stick-fixed neutral point = center-of-gravity location at which 
06 ¢/OV = 0 in steady one-g flight 

Stick-fixed, maneuvering neutral point center-of-gravity 
location at which 06 e/On = 0 for a constant-speed maneuver 

Stick-free neutral point = center-of gravity location at which 
O( F/qg)/OV = O in steady one-g flight 

Stick-free, maneuvering neutral point = center-of-gravity loca 
tion at which OF /On = 0 fora constant-speed maneuver 

Aerodynamic center = center-of-gravity location at which 


OC\,/0a@ = 0, with speed being held constant 


INTRODUCTION 


6 ie HANDLING CHARACTERISTICS OF AIRCRAFT were 
subjects of extensive study even before man’s 
first flight. From the earliest days of aviation, the 
designer's central problem has been to develop work 
able controls and to supply stability or, at least, to 
hold instabilities within controllable limits. As the 
art of designing and building aircraft progressed, pilots 
became more discriminating and demanded airplanes 
with better flying qualities. Eventually, the desirable 
handling characteristics became well enough defined 
to warrant writing specifications outlining minimum 
requirements. Some of the most important of these 
requirements relate to longitudinal stability. Since 
static longitudinal stability decreases as the center of 
gravity is moved aft, it follows that the longitudinal 
stability requirements can limit the aft c.g. location. 
However, inappropriate requirements may be detri 
mental to the utility of an airplane if they impose 
unnecessarily restricted loading limits. 

Before the advent of jet aircraft, the requirements 
for satisfactory longitudinal flying qualities were well 
established and generally accepted. These same re- 
quirements—couched in terms of static margin, rela- 
tive to the neutral point-c.g. relationship—are em 
ployed today. The curve of elevator-to-trim versus 
speed-in-flight is presently used to evaluate the sta- 
bility characteristics. Positive stability is indicated 
when an increment in down elevator is required to trim 
at increased speeds. However, compressibility and 
aeroelastic effects may introduce reversals in the trim 
curve with speed—just as power effects do with a pro- 
peller airplane—even though the airplane is very stable 
to disturbances in the angle of attack at a fixed speed. 
Nevertheless, the airplane may be quite acceptable for 
the pilot to fly under these conditions. Since this is 
true, it appears timely to reconsider the basis of the 
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” ‘ie AXIS 


Fic. 1. Coordinate system 


longitudinal stability requirements with an eye toward 
developing more fitting criteria. 
THEORETICAL CONSIDERATIONS 


The longitudinal motion considered will be restricted 
to the plane of symmetry and will allow freedom in the 
angle of attack, pitch angle, and speed, as shown in 
Fig. 1. The wind axis system is employed. The 
elevator will be assumed to be fixed for this part of 
the analysis. 

The stick-fixed equations of motion become: 


Aa(T, — dD.) i AT 7) - Dy) —_ 
(W /qS) (A@ — Aa) = (Wg) [V/(qS)] 


AaL, + AVLyr = (W/g) [V (gS) | (0 — a) 
AaM, + 4M, + AVMy + 6M, = [I,/(qS) 0 
where 


Do = {I (qS) | (OD 0a); My = 
[1 (gS.)] (OM 06); ete. 


This system of equations may be expanded to yield 


the characteristic equation of the system: 
Ad! + Br? + Ch? + DA+ E = 0 (2) 
where 


A = (C_/g)*V[J,/(q5-) | 


B= —(C, g)?V(Ma + My) + 
(C. g) Ul, (gS) [La -— ViTy -— Dr) 
C = —[I,/(q¢Sz)] [((Da — Cr)Ly + 


(Ty; — Dy)La] + (Ty — Dy) Ci(V/g) (My + Ma) - 
(C,./g) [LeMy + Ci(V/g) Me] 





D = (Ty — Dy) [L-Ms + Ci(V/g)Mo] + 
Cr [My-Ch( V g) = LyM, | -F 
(D. — Cr) [LvMg + Ci(V/g) My] 


k= CLL My _ Ly M,) 


Neutral static stability is indicated when - = 0. 
A positive value of / ensures that a pure divergence 
does not exist, while a negative value of / permits a 
positive real root, indicative of a pure divergence. 

It is evident that the change in moment with speed 
(My), as well as the change in moment with angle of 
attack (/,), influences the static stability. A negative 
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TABLE | 
Airplane Characteristics 
Airplane 1 Airplane 2 
Characteristics Propeller Jet 
V 89 knots 528 knots 
] 26.6 lbs. per sq-ft 945 
7 —(0). 0020 0 
"D 0.0936 0.0250 
dD, 0.00125 QO OQOOOTOS 
Dea 0.634 0 
0.0124 QO. O00125 
me 0. 0556 
C; a 5.43 $ 59 
Wy, —0.00125 —(0, OOO03B2 
T.’ (Total 0.30 0.0250 
c.g 31 per cent MAC 44 percent MAC 
Ma —().309 —().115 
Va —0.629 —(). 0514 
M x —0.346 — (0.00922 
GW $2 200 Ibs L100 ,000 Ibs 
I; 500 , 000 slug-ft.* 1.65 X 10 slug-ft? 
> 1,408 sq ft 1,900 sqft 
( 13.08 ft 20.7 ft 


value of ./, and a positive value of .1/;- are in the stable 
sense. For stable conditions then, increased angle 
of attack and increased speed result in diving and 
stalling moments, respectively. 

The quartic equation governing the pitch motion 
generally yields two oscillatory modes that make up 
the airplane’s response to a disturbance. The short 
period motion is generally heavily damped, while the 
long-period motion—classically called the phugoid 
motion—is ordinarily lightly damped. In the classical 
phugoid motion for which ./;- = 0, the phugoid period 
is a linear function of speed, theoretically given as 
P = 0.1381;,;.. The short-period mode is then pri- 
marily a function of 1/7, and is usually little affected by 
My, except at center-of-gravity locations which have 
AM, values near zero. 

If these two motions are stable and adequately 
damped, the airplane must be considered to have satis- 
factory stick-fixed longitudinal stability. The sta 
bility requirements of reference 3, however, are not 
expressed in terms of the damping of these motions 
The longitudinal stick-fixed stability requirements are 
expressed therein in terms of the stick-fixed neutral 
point. The stick-fixed neutral point is defined as the 
c.g. location at which the elevator position gradient 
(with changing speed in a steady one-g flight with fixed 
throttle) is zero. 

Mathematically, this condition requires dC,, dCz in 
steady flight to be zero, but 


dCn/AC, = [(OCn/Oa) (OC, Va) ] + 
(OCy Ol” 


Ol’ 0C;) (3 
In steady one-g flight, 01° OC; = —I (2C;), so at the 
neutral point, 


[(OC, Oa) (OC; Oa) | — 
(OC, OV) [V (2C 
It should be noted that in maneuvering flight at con- 
stant speed, 0V’/OC, = 0, and the effect of a change 
in pitching moment with speed cannot influence the 
total value of dCy/dC,. It is then possible to have 
the elevator deflection per g in maneuvering flight 
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CRITERIA FOR LONGITUDINAL 


at center-of-gravity locations well aft of the 
Ol has a large negative or un 


stable 
neutral point, if OCy, 
stable \ ilue. 

The dynamic stability with free controls is also of 
sreat importance. The stick-free longitudinal motion 
involves four degrees of freedom with the elevator mo 
additional degree of freedom 


tion prov iding the 


tbove those considered thus far. For a mass balanced 
control system, the hinge moment equation is given 


below 
Mall, 61, + Asi; ++ 6H; = 


J, (gS.C, 


all, + 
} 6+ 6 oe 


Equations (1) and (5) will then define the airplane 


dynamics with the stick-free. 


ANALYSIS OF REPRESENTATIVE AIRPLANES 


To clarify some of the principles involved, two repre- 
sentative airplanes will be considered. The first is a 
propeller-driven airplane flying at S89 knots in a power- 
approach condition, the second, a jet airplane flying 
at Mach 0.8 in the tuck range. The airplane char- 
acteristics are tabulated in Table 1. 

The stick-fixed, stability quartic A\' + Bd* + 
C\? + DX + E = 0 may be factored into two quadratics 
e+ A, + B,) and (A? + CA 4+ D,). 
ratics are representative of the short- and long-period 


These quad- 


If complex roots denoting an oscil- 
i.e., A2 = a + 1b—the 


stick-fixed modes. 
latory mode are assumed 
period and damping of the mode may be determined 


from the values of a and bd. 


T.,, = 0.693/a 
and 
P = 6.28/) 


It is easily seen that the cycles to one-half ampli- 


tude may be determined as 
C.,, = 0.1104 b/a 


The quadratic factors are functions of the center- 
of-gravity position, and, for certain locations of the 
c.g., the roots of the quadratic may both be real num- 
bers. The motion is then nonoscillatory in nature and 
isa pure divergence if one or more roots are positive. 

The roots of the quadratic equation defining the 
stick-fixed short-period motion for both airplanes are 
plotted on the complex plane in Fig. 2. Inspection of 
this figure shows that, as the c.g. is moved aft, the 
imaginary values of the roots disappear, and one real 
root eventually moves to the right-hand half of the 
plane where a pure divergent motion is indicated. A 
pair of lines corresponding to the damping to one-half 
amplitude in one cycle are also included. It may be 
seen, as the center of gravity is moved forward, that 
the roots move farther from the real axis and eventually 
will cross the C,,, = 1 lines. Outside these lines, the 
damping of the short-period oscillations is poorer, due 
to excessive static stability. Past experience indicates 
that the short-period pitch oscillation should damp to 
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one-half amplitude in one cycle or less. Therefore, 
certain limits to the roots of the short-period mode are 
specified. 

However, a further stability requirement for the 
short-period mode is needed. The return to the initial 
condition after a disturbance should be reasonably 
rapid. As a tentative value, based on present experi 
ence, it is proposed that the time required to damp to 
This 
requirement will prevent an excessively sluggish air 
On the complex plane, this excludes 


one-half amplitude should not exceed 1.5 sec. 


plane response. 
a real root to the right of — 0.46. 

The two requirements for short-period stick-fixed 
stability—namely, (1) damping to one-half amplitude 
in one cycle or less and (2) damping to one-half ampli 
tude in 1.5 sec. or less—will then define an area on the 
complex plane for satisfactory values of the roots of the 
stability equations. 

It is apparent that these requirements will define the 
aft limit for the center of gravity of the airplane. For 
the propeller-driven aircraft under consideration, the 
For the 


aft limit is indicated at 32 per cent MAC. 
jet airplane flying at Mach 0.S at sea level, the aft 
limit is at 48.2 per cent MAC. 

A comparison of the aft c.g. location with the aero 
dynamic center, c.g. for divergence, stick-fixed maneuv 
ering neutral point, and the level-flight neutral-point 
is shown in Table 2. 

It is significant to note that, for Airplane 1, this 
proposed aft c.g. location is 4.4 per cent MAC forward 
of the aerodynamic center or c.g. at which OC, 0a = 0, 
while on Airplane 2, the dynamics criteria aft c.g. is 

TABLE 2 
Comparison of Criteria for Aft Center-of-Gravity Location 


Dynamics 


Aerodynamic Stick-fixed criteria Diver- 
center maneuvering aftc.g gence Neutral 
Airplane WC yw Va 0 neutral point stick-fixed c.g point 
1 (Propeller 36.4 48.4 32.0 35.0 28 1 
2 ( Jet 46.5 95.0 18.2 $9.5 20.8 


All c.g. locations are per cent M AC 
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rr TW ti LITT TL ae ations, a positive Cy, can appreciably improve thy stl 
L++4+—_+-+4, Ha aft center-of-gravity limit, as shown by the calculations sir 
ettiteo or,” a for Cua = 0.0020. wit] 
BSaesee NUMBERS INDICATE Consider the stick force-per-g in a steady pull-up for 
Sin ae es C.G LOCATION : ; : é ; airp 
a6 a . the two cases with Airplane |. For this aircraft the 
ee cl TTT minimum acceptable stick force-per-g is 22.5 Ibs.-per-g T 
> | ro 1 #8 The aft-most c.g. at which this stick force-per-g re bail 
Zo 138 quirement is met is 33.1 per cent MAC for Cy, = () T 
3 1H | rT and 45.5 per cent MAC for Cy, = 0.0020. For the mot 
- +++ case in which Cy, = 0, the stick-free damping require asl 
CTT {| | ments limit the usable aft center-of-gravity location ios 
T Hi Co But, with Cy, = +0.0020, the stick force-per-g and sis 
-8 | ae we stick-free dynamics are satisfactory to approximately airp 
‘SATISFACTORY | 15 per cent MAC. This is appreciably aft of the stick a 
£0 cq ci fixed divergence limit of 35 per cent MAC. oscil 
In the author's opinion, center-of-gravity positions mod 
7“ — - - aft of the stick-fixed divergence point are unsafe, even _— 
woe ; though the stick-free dynamics and stick force-per-g the | 
Pes. 3. Rovts for stick-teee short-perted motion are satisfactory. An extremely dangerous condition cenit 
would arise if, with this c.g. position, the aircraft feel How 
located at 1.7 per cent aft of the aerodynamic center. is lost somehow, as when there is binding of the con- be s 
Flying at 89 knots, Airplane | has a stable dead-beat trols or failure of the synthetic feel-system. Conse long: 
short-period oscillation over a c.g. range from approxi quently, the safeguard of a dynamically stable aircrafi fied 
mately 28 to 35 per cent MAC but has very sluggish in the stick-fixed condition is deemed necessary } later. 
response near 35 per cent MAC. Airplane 2, flying With this one restriction, however, it is considered | alogy 
at Mach 0.8, has a stable dead-beat oscillation over that the stick-free dynamics are of more fundamental | the « 
a c.g. range from approximately 45.5 per cent MAC importance. vet e 
to 49.5 per cent MAC. Note also that  stick-fixed The premises upon which the aft center-of-gravity “Tt 
divergence occurs well forward of the stick-fixed man- limit is set should be based primarily upon stick-free longi 
euvering neutral point. This is due to the unstable considerations. It is believed that control force, ta 
values of O0Cy OV for the airplanes used as examples. rather than control position, sensory cues are of over- eis 
It is evident that the neutral point has little correla- riding importance in the pilot’s opinion of aircraft denies 
tion with the airplane short-period motion for these air- flying qualities. In summary, from consideration of ticule 
planes with their unstable OCy/OV characteristics. both stick-fixed and stick-free flight, the following re- | to ro 
At the neutral point, both airplanes have very strong quirements for short-period longitudinal stability are tione 
short-period restoring tendencies. Airplane 2, in par- proposed. The aft center-of-gravity limit: (1) shall seque 
ticular, is extremely stable from a short-period stand- not be behind the c.g. location at which the stick force- is los’ 
point, the aerodynamic center being 25.7 per cent per-g is less than a specified value, depending upon | same 
MAC aft of the neutral point. It is clear that the neu- airplane type and limit-load factor; and (2) shall be | the fc 
tral point is not a reliable means for determining the forward of the stick-fixed short-period divergence poimt. | jog 
aft center-of-gravity limit, and a specification that sets The short-period, dynamic stability characteristics | he 
the aft c.g. limit as a function of the neutral point may | long- 
be unnecessarily conservative. Conversely, however, We SRERER RSE EE DEER BERT by th 
for aircraft with stable values of O0C\,/OV’, the neutral py tt fp tt dition 
point could be an unsafe criterion for the aft center-of- CI t ey eS eS Oe we ] eee e 9 eas Howe 
gravity limit. a 28. “SE SSSS 58 25 shoul 
In any discussion of the short-period, longitudinal a <= : tt + condi 
dynamics, it is, of course, necessary to consider the = TT T oe ia gence 
stick-free case. Calculations of the stick-free dynamics | ew tude i 
and stick force-per-g have been made for Airplane 1, 2 a In’ 
and the results are shown in Figs. 3 and 4. Two com- rd | tt of th 
binations of hinge-moment characteristics were used: 5 weet Comp 
Cua = 0 and Cy; = —0.0040 per degree; also Cy, nh Zzuse long-4 
+().0020 and Cy; = —0.0040. For the case in which me and < 
Cua = O, the center-of-gravity location at which 7), = SH SRASSESR Ae SS PARR SS SE” center 
1.5 sec. is 31.8 per cent MAC; when Cy, = 0.0020, the a | or a ie phugo 
comparable c.g. position is approximately 50 per cent SLES ESATA eS feipial Lo becor 
MAC. It is seen that for the zero-floating-tendency — Oo 3 “40 50 60 quite 
condition, the stick-fixed and stick-free dynamics are CENTER OF GRAVITY-% MAC ployec 


nearly identical. However, from dynamics consider- Fic. 4. Characteristics of stick force-per-z, Airplane | range 
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CRITERIA FOR LONGITUDINAL 


shall be such that: (1) the stick-free, short-period, 
airplane oscillation damps_ to one-half amplitude 
within one cycle; and (2) the stick-free, short-period 
airplane oscillation damps to one-half amplitude within 
] sec. 

These latter requirements may also limit the aft 
loading. 

The characteristics of the stick-fixed, long-period 
motions are not closely related to center-of-gravity 
location for either of the airplanes considered. Table 3 
shows the characteristics of this mode. The main 
characteristics in common between the two example 
airplanes are that, as the c.g. is shifted aft to produce 
an unstable short-period motion, a stable long-period 
oscillation comes into existence; and the long period 
mode is most unstable, in both cases, in the area just 
forward of short-period divergence. It is apparent that 
the long-period mode cannot be used as a condition to 
limit the location of the airplane’s center of gravity. 
However, it is believed that some requirement should 
be specified, relative to the rate of divergence of the 
long-period mode and in a manner similar to that speci- 
fed for the rate of divergence of the spiral mode in 
lateral-directional stability considerations. The an- 
alogy between the long-period, longitudinal mode and 
the spiral mode is strong. Both can be unstable, and 
vet easily controlled, if the rate of divergence is low. 

It is entirely possible that the stick-fixed, long-period 
longitudinal requirements may be eliminated in favor 
of stick-free equivalents. The pilot flies almost ex- 
clusively by sensations of force with respect to speed 
changes. To slow an airplane down from any par- 
ticular trim speed, the stick must first be pulled back 
to rotate the airplane nose upward, and then reposi- 
tioned for stabilized trim at the desired speed. In this 
sequence of operations, the original position reference 
is lost, and only the force trim reference remains. The 
same reasoning holds true if a speed build-up is desired; 
the force reference remains and the position reference is 
lost. 

As previously mentioned, since a mildly divergent, 
long-period, pitching oscillation can be easily controlled 
by the pilot, it should be permissible under such con- 
ditions as take-off, wave-off, and high-speed runs. 
However, a limit to the allowable divergence rate 
should be set. It is recommended that, for the flight 
conditions under which a stick-free, long-period diver- 
gence is tolerable, the motion shall not double in ampli- 
tude in less than 8 sec. 

In Table 4 are presented the results of the calculations 
of the long-period, stick-free motion for Airplane 1. 
Comparing the results for Cy, = 0 with the stick-fixed 
long-period dynamics, it may be seen that the period 
and damping characteristics are very similar. Aft 
center-of-gravity movement first results in decreased 
phugoid stability, and then, as the short-period motion 
becomes unstable, the long-period motion becomes 
quite stable. When a value of Cy, = 0.0020 is em- 
ployed, the long-period motion is stable over the c.g. 
range of the calculations, although it becomes just 
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mildly convergent at a c.g. location of 41 per cent 
MAC. The very favorable effect of the positive Cy, 
is apparent. For both cases of Cy, = O and Cy, 
0.0020, the stick-free phugoid divergence-time exceeds 
the proposed time to double amplitude——of 8 sec. over 
the c.g. range—which is permissible from short-period 
considerations. Thus it is satisfactory within the 
framework of the proposed requirements. 

However, in the cruising and climbing conditions, it 
is necessary that the airplane not require constant 
attention from the pilot to maintain the desired flight 
attitude. This requires positive stability under these 
conditions. Consequently, it is believed that the 
requirements should demand a pull force to stabilize 
at a decreased speed and a push force to stabilize at 
an increased speed in cruising and climbing flight. 
This force gradient should be of sufficient magnitude 
to maintain essentially the trim speed after a sharp 
control displacement and release. 


DISCUSSION 


In essence, it is proposed that the characteristics of 
the airplane control forces and short-period motion, 
rather than the neutral point, be used as the criterion for 
the determination of the aft loading-limit. 

The two airplanes considered have unstable vari- 
ations of pitching moment with speed, and these exert 
a large influence on the neutral-point location. In 
Airplane |, this is due to power effects; in Airplane 2, the 
cause is compressibility. Because of the unstable 
OCy/Ol, the neutral point is very far forward of the 
center-of-gravity locations where the restoring tendency 
of the airplane to a pitch disturbance is small. 


TABLE 3 
Long-Period Pitching Motion—Stick-Fixed 


Center of gravity, 


Airplane per cent of MAC T1/; Period r; 
l 21 15.4 58.2 
26 62.0 
31 8.10 
36 2.01 29.3 
41 2.52 21.2 
9 34 34.0 221 
44 33.6 225 
47 27.3 95.0 
50 8.1 41.3 
54 29.5 74.7 
TABLE 4 


Long-Period Pitching Motion—Stick-Free 


Center of gravity, 


Airplane per cent of MAC T\/; Period T; 
Cua = 0 
l 21 15.9 56.1 
296 63 
31 8.) 
36 2.06 26 .2 
4] 2.62 20.5 
Cua = 0.0020 
91 18.6 35.8 
26 18.0 39.2 
31 16.9 15.6 
36 15.3 62.5 
41 155 
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The long-period mode, which may become unstable 
due to adverse variations of pitching moment with 
speed, does not necessarily become more stable as the 
c.g. is shifted forward. Therefore, it does not appear 
to be an appropriate criterion for the determination 
of the aft center-of-gravity limit. However, the 
motion can be annoying and unacceptable if its rate 
of divergence is high. Because of this fact, it is con- 
of the long-period pitching motion. It is furthermore 
considered appropriate to require stick-free, long- 
period, longitudinal stability in cruising and climbing 
flight. 

The fact that the short-period motion is of overriding 
importance to an airplane’s flying qualities has also 
been demonstrated in the free-flight model tests re- 
ported by Schade.' In these tests, center-of-gravity 
shifts, and thrust-line offset were used to obtain vari- 
ous neutral-point locations for particular values of the 
aerodynamic center. The model flight behavior was 
then evaluated for numerous combinations. In Fig. 5 
the results are re-plotted from reference 1. 

It is readily seen that the neutral-point location 
had negligible effect upon the flight behavior, as rated 
by the pilots. The important factor was adequate 
constant-speed stability to produce the strong restoring 
tendency in the short-period motion. Considerable 
negative static margin—1in terms of neutral point, — 11 
per cent MAC-—still permitted good flying qualitites 
if enough aerodynamic center margin was present. 
This aerodynamic center margin is measured by the 
value of (OCy/0a)/ (OC, 0a). 

The requirements that have been advanced are con- 
sidered to be sufficiently stringent and comprehensive 
to ensure an aircraft that is safe and one that should 
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be readily accepted by pilots. However, an attempt 
has been made to formulate the requirements in 4 
manner that stipulates no unnecessarily severe or jn- 
appropriate conditions. Although no claim is made 
that these proposals are inflexible, it is believed that 
the concept of the specification requirements which 
has been advanced is important and that an approach 
of this type will provide better aircraft. 


CONCLUDING REMARKS 


An attempt has been made to improve the longi- 
tudinal stability requirements. This problem is con- 
sidered of importance, since a poorly founded specifi- 
cation requirement can severely limit the utility of an 
aircraft. Specifically, it has been shown that the 
steady-flight neutral point has little significance in de- 
termining the allowable aft center-of-gravity location. 
The short-period, dynamic characteristics and the stick 
force-per-g are considered to be the important factors 
in determining the allowable aft c.g. location. 

The recommended requirements are summarized 
below. The aft center-of-gravity limit: (1) shall not 
be behind the c.g. location at which the stick force-per-¢ 
is less than a specified value, depending upon airplane 
type and limit-load factor; and (2) shall be forward of 
the stick-fixed, short-period divergence point. The 
short-period dynamic stability characteristics shall be 
such that: (1) the stick-free, short-period airplane 
oscillation damps to one-half amplitude within one 
cycle; and (2) the stick-free, short-period airplane os- 
cillation damps to one-half amplitude within 1! » sec. 

These latter requirements may also limit the aft 
loading. 

The long-period, stick-free dynamics may be un- 
stable, except in cruising and climbing flight, but the 
motion shall not double in amplitude in less than § 
sec. In cruising and climbing flight, positive-force 
stability must exist so that push and pull forces are 
required to stabilize at speeds above and below trim, 
respectively. This force gradient should be of suffi- 
cient magnitude to maintain essentially the trim speed 
after a sharp control displacement and release. 

Some of the proposed requirements have already 
been stated in the Air Force and Navy Handling 
Qualities Specifications* and thus conform to generally 
accepted practice. These specifications, recently un- 
classified, have provided the basis for flying qualities 
design of military aircraft from 1948 through 1954. 
It is clear that the newest concepts introduced in the 
present pages must be confirmed by specific flight tests. 
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The Effect of a Finite Hub on the 
Optimum Propeller’ 


B. W. McCORMICKt 


The Pennsylvania State Unwersity 


ABSTRACT 
Goldstein's analysis of a propeller satisfying the Betz condi 
tion is extended to the case where the hub radius is appreciable 
relative to the propeller radius. The solution, which satisfies 
the boundary conditions of an infinitely long cylinder representing 
the hub, is expressed as an infinite series of modified Bessel func 
tions It is shown that to a close approximation the differences 
between the present case and that with a zero hub radius can be 
calculated without too much difficulty. Examples are pre 
sented which show that for a trailing vortex sheet of given pitch, 
the insertion of a hub increases the value of the bound circulation 
of the line vortex generating the sheet This change is more 


pronounced the smaller the number of blades 


SYMBOLS 


= coefficients in the I'-distribution for 7, Zr ZR 
A = coefficients proportional to a,,’s 
= coefficients in the ['’-distribution for 7, Zr ZR 


number of blades 


( = coefficients in the I’-distribution for r > R 
fm = function defined by Ea. (14 

Sn ( oe = function defined by Eq. (18) 

D(a = modified Bessel function of the first kind 

K(x modified Bessel function of the second kind 
= index 

, = index 
= radial distance from propeller axis 

R = propeller radius 

> = Lommel’s function 

Mm = function defined by Eq. (21) 

T\. ,(@) = function defined by Eq. (22) 
= advance velocity of the propeller 
= advance velocity of the helical vortex surfaces 
= axial distance behind propeller 

I = bound circulation 

€ = small correction to 4,, to obtain exact value of a,, 
= independent variable defined by Eq. (4) 

6 = angular coordinate 

re = independent variable defined by Eq. (4) 

¢ = velocity potential 

oi) = portion of ¢ as defined by Eq. (11) 

w = propeller rotational velocity 


Subscripts 
= 2 


() = pertaining to quantities defined at r = R 
h = pertaining to quantities defined at the hub radius 
Superscript 

Note 
of the function is given by the product of the unbarred 


a bar over an argument indicates that the argument of 


argument with the order of the function 
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INTRODUCTION 


i bw PROBLEM TO BE CONSIDERED is that of calcu 
lating the distribution of bound circulation of an 
optimum propeller, the hub radius of which is rela 
tively large compared to the propeller radius. By 
optimum propeller is meant a propeller in potential 
for minimum 


flow which satisfies the Betz condition 


energy loss namely, that the generated trailing vortex 
sheet, in the steady case, moves as a rigid helicoid. 
The limiting case of this problem is that of the optimum 
propeller with a zero hub radius. The classic solution 
by Goldstein to this problem can be found in reference 
1. The present development represents an extension 
of Goldstein’s methods to the more general case of a 
finite hub. It will be shown that the general case is 
most easily calculated by considering the difference 
the 


which are available for the zero hub case can be applied 


between it and Goldstein's case. Thus results 


directly to the more general one as developed here. 


Formulation of the Potential Problem 


The formulation of the potential problem is nearly 
identical to Goldstein's case with the exception of the 
additional boundary condition at the hub surface. If 
R is the propeller radius, w the angular velocity, and 
v the velocity of advance of the propeller, the equations 
of the trailing vortex sheets are 

0 — (w/v) = 0, (270 /B (= 13.3.0. (10 


forr, Zr Z R, where 


r, 6,2 = polar cylindrical coordinates 
B = number of blades 
1 = hub radius 


Eq. (1) assumes the propeller to be lightly loaded with 
the effect of the induced velocities upon the position 
of the trailing vortex sheets being neglected. 

If w is the velocity of advance of the trailing helical 
vortex surfaces in the z-direction and if ¢ is the velocity 
potential, then the boundary condition that the ve- 
locity normal to the surface of the sheets be zero is 


wre = wr(Od Os) — (vr) (Od O88 (2) 


for 


6 — LE 


(ws/v) = 0, (27n/B) % £ 


At the hub surface the condition is imposed that the 
normal velocity must vanish. If the hub is replaced 


by an infinitely long circular cylinder of radius r,, 
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the imposed boundary condition means that for all 6, 


it must hold that 
(0¢/Or),-,, = O (3) 


In addition to conditions (2) and (3), grad @ must vanish 
as ry becomes infinitely large. 

From its symmetry, the problem can be reduced to a 
two-dimensional one with the independent variables 


being 


C = 6 — (we/v 
™ (4) 


In terms of these new variables, the boundary condi 
tions can be rewritten as 


(Op OF) 2en = —[u?/(1 + 2*)] m Zr ZR (5) 
0 


(Od Op), = 9 (6) 


where yp, = wr),/v. 

These boundary conditions take the constant quantity 
w/w temporarily to be unity. The actual ¢ will be 
given by the product of the ¢ thus obtained and the 
quantity wv /w. 

The differential equation to be satisfied by ¢ is 


V*@ = 0 (7) 
which in terms of u and ¢ can be written as 
[u(O/Ou) |?*@ + (1 + uw?) (0°f/0¢?) = O (8) 
where the operator 
[u(O/Op) |? = ul(O Ox) [uC Oz) } 


In addition to satisfying Eqs. (5), (6) and (8), the 
function ¢ must be single valued, have vanishing deriv- 
atives as uw becomes infinitely large, and be continuous 
everywhere except at the vortex sheet. 


Solution forr>R 
For radii larger than the radius of the propeller, 
the solution of ¢ will take the same form as Goldstein’s 


case—namely, 


2C, Ke,(Bnpz) 


: sin Bn ¢ 
n=1 B Kp,(Bnw) 


40) = 
(9) 


In the above equation, Az, (Bn yu) is the modified Bessel 
function of the second kind and the C,’s are undeter- 


mined constants. The constant yo is defined by 


Mo = wR/v 


At this time a matter of notation is emphasized. In 
order to avoid unnecessary writing, a bar notation will 
be used for the argument of any Bessel function to 
signify that the argument of the function is given by 
the product of the term under the bar with the order of 


the function. For example, 


Tia) = I,(vp) 
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Using this bar notation, Eq. (9) will become 


26. Kp,,(B) . 
= > : sin Bunt n ae ae! 
n=1 B Kx,,(io) 


10 


The constants, C,, will be determined later from th: 


conditions that ¢ and its derivative be continuous 4} 


Kh = Mo. 


Solution forr, <r <R 

In this region the solution for ¢ again follows quit 
closely the development presented by Goldstein. Hoy 
ever the differences between the two are sufficient t 
warrant a detailed presentation. 

Following the lead of reference |, the velocity poten 


tial is assumed to be of the form 


From Eq. (8), it follows that ¢ must satisfy the follow 


ing differential equation: 


[u(O Ou) }?o1 + (1 + w*) (0%, 06?) = 
[u(O/Ou) |? [w?/(1 + w]e (2 


Since @ is periodic with ¢ only the region for which 
0 Z ¢ Z(27 B) need be considered. In this region, 
¢ can be expanded in a half-range cosine series 
. S ys cos [(2m + 1)/2|Be , 
= — , i. » 
B 7B m=0 (2m + 1)° 


¢; is expanded in a semi-infinite cosine series in ¢ for 
this interval 


2m + | 


od: = fo(u) + z, fm(u) COS Be l4 
) 


m ( - 
For Eq. (12) to be satisfied for all values of ¢, it follows 


that 


[u(O/Op) |?fo(u) = (7 B) [u(O Ow) ]? [w?/(1 + w*)] (15 
and 
se [Im +1 _|? 
( Fcdl )— (1 a =) B| fm (um) 
r 7 oe ie ae ' 


S | Oo \- pe 
L 16 
7B (2m + 1)? Ou 1 + pe 


Eq. (15) can be integrated at once to give 
fo = (w/B) [w?/1 + w)] + Ging 17 


where C, is the constant of integration. The second 
constant of integration was taken to be zero since the 
level of ¢ is arbitrary. 

A new function g,,(u) is now defined by 


| 2 7 
S ! “= | (18 


Feat ; = ; . LZm\b 
ele rB (2m + 1)? 1+ pw 


In terms of g,,(u) Eq. (16) becomes 





then 


A pal 
funct 
solut 


the s¢ 


wher 


term 


? sin 


rhe 


of E 


pelle 


1M 


ad ~~ p- p-( p — » 
: ) g,(zZ + T(z = |]— oe = 
az M he 
{ 2m + | Ze v*(y? — 27) (yp? — 4+ 
| 2° — 2 B  m(Z) s* (iS = nae ve ~ K, (i 23 
L - j Ta 
where 2 i}(2m + 1) 2]Bu. The general solution for g,,(u) includes the homo 
geneous solution of Eq. (20) which is a linear combin 


f V, is an operator signifying . ‘cia : : 
om th ” iliac ation of the modified Bessel functions of order [(2m 
it P - 
10US at v, [s(d ds) |? + 2° ‘ 1) 2|]B and argument g. Thus g,,(u) is given by 
1 Ea. (19) re: 2 mi TL omti . (2) + Amlom 
then Eq. (19) reads £m (Me 1, 2m+1 , B 2m+l 
- o ») 7 , 
Vez ' p) £m \2 . (20 er Z 24 


(2 


ei 





S quite 
How A particular solution of the above equation is Lommel’s Redefining the coefficients a,,, b,,, and Cy to include the 
lent to function discussed in reference 2. If this particular values of the modified Bessel functions for argument 
solution is denoted by S_ 5,,,4.; (z) then one form which of @ and the number of blades, expression (11) for @ 
1 ; 
poten ‘ becomes 
the solution takes 1s 
m 1 | : 8 T am+i, Gi 2a lom+i , (i 
l] Br o » ss 2 aA 2 
2 m=0| 7B (2m + 1 B fois vi 
S 2m4+1 , 2 xX - B 
ollow B cs 2n + 1 7B 5 
~ sin 2 9 A om+1 . (ii) 
= e 2b, 5a 8B 2m + 1 iv . 
/ (za f fas cos : Be + CG 1n (25 
= ' B . i‘. = 'B P BiKoy, 4) B \M - A 
) = 2 
2 2m+1 xB 
2m l ( 5 ) ; " e ‘ . — 
Be - Kom 3 y (21 In the case of a zero hub radius, the coefficients, ),,, 
° ) ; 
which sz are identically zero since the A»,,,; — (@) functions are 
egion, : 3 —~S 
where ata , , 
| infinite for zero argument. In the general case of a 
x x* finite hub however, these terms are retained. 
‘ 92 — y’ (22 — y?) (42 — yp? The a,,, 6,,, and Cy in the above equation and the 
” 
6 C, in Eq. (10) are to be determined by the three con 
\ | ; 
9 79 , 22s ditions: 
<e (2 - pn) ( — ({ 
© for 2 V | v ) l 
' bladed , ' ie (1) @ continuous atr = R 
‘or two, six, ten . bladed propellers the coethcient . 
For t _ é ‘ “e I ; ' 7 (2) grad @ continuous atr = R 
> Kon 4 (a) is purely imaginary and can be 
4 of the 2m- i. LL pu i id 7 (3) O@ Ou = (Qatr=yr, 
dropped. Goldstein defines Eq. (21) without this 5 . ; i E 
PI | , 
: Fe Evaluation of Undetermined Constants 
lows temas T 5,4 1 (ft). 
. l nd B - . . 
2 In order to evaluate the coefficients a,,, 5, Co, and 
C,, it is first necessary to express @ for r < R in terms 
as) | 7 3 S 2m+1 , A) — of the sine of multiples of ¢. This is accomplished by 
ee ai expanding cos |[(2m + 1) 2]|B¢ in a half-range sine 
) 1 - ; . .  * : ; ‘ 5 
<m | 2 2 ) K (a) = series. If this is done then it will be found that in 
2 =—° = order to satisfy the condition that @ be continuous 
(2n |) 2|Bx atr = R for all ¢, the following must hold. 
1. (@) —T? _. ) C2 
. [2m + 1 Br I = B 1 se 
2 sin om+1. (h 
16 ) 5) SoG S I 2a b 
a T aie x 
—— ; : ) ; T 0| tB mm + | B B 
he function 7 .,,4) — (@) is also a particular solution 
1, -—— B n 2¢ 
— ; 7 ; 26 
. ot Eq. (20). For four, eight, twelve . . . bladed pre (2n (2m + 1 B 
l ele wets 
: pellers, the form (21) for S_ 5,,4, . (@) is indeterminate. ee — 7 
“le Similarly, to satisfy the condition that grad @ be con 
| a c z . bs ‘ 
a . . ° ° . . ° . 
2 In this instance the asymptotic solution for .S) (ax tinuous at ry = Ra second relationship is obtained be 
tne . . . . = 
, given in reference 2 must be used. tween the coefficients. 
, ~" , oi , , = 
18 TD 2m+1 , (io) Tom+1 , op amt, Uh ' 
v 1(2m + 1 S ie 2a, B 2b. 3 C, Ap,’ (i ” 
= T _ T 2a = uy 
0 4n (2m + 1 cB (2m + 1) B fons, . @ BD hs. 1, B Kp, (i 
= B ; 
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A prime in Eq. (27) denotes the derivative of the function with respect to the argument g@ and not the independent 


variable u. 


Finally the condition that 0¢@ Ou = 0 at the surface of the hub results in the following relationship between ,, 
and b,,. 
Kom+1. (io) Toms.’ LP oom4i , (ih 
+ B B | = B 
b. = . : a 7. : 29 
K 1, An) Coe (Po) T (2m + 1 
Eliminating b,, and C,, between Eqs. (26), (27), and (29) results in an infinite number of equations (7 = 1, 2, 3 
for the infinite number of unknowns a,, (m = 0, 1,2...). 
x Lom +1 i (fo) Kom +1 ; (Ho) l i (Ey } 
n in 5 B > E : B 
= = Veme = 1b) — ; = ; — 
4 m=0 (2n)? — (2m + 1)? | Tom41 . (io) Roms.’ (in) Lom41 (i) 
: B : B : B 
=, K», l , (Ho) Dom 41 i (En) , +, T 2m+1 (fy) 
Kp,’ (io) Watliad > 5 Kp,’ (io) oe 
Bn \b 2 2 Bn \b 
22 — i= — > = 2n — pm oS Xx 
Kp, (io) (om 1. (An) Bon, bs (fio) | Kpn(fo) m=0 (2m + 1)*[(22)? — (2m + 1)?] 
> > 
Kom ot (Ko) aL 2m+ (y,) ij 2m+i : (fo) 
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ae ee meee Gat oe oer o— —e 
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7 > 
Kk . 1 e ( Ho) 7 Dy, 1 E ( My 
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Ko, 1 (fy) 1 wad ( 
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The effect of the finite hub is evident in the above 
system of equations. The terms in the square brackets 
represent the difference between the general case and 
the zero hub case. Each of these brackets will be 
close to unity. This will be especially true, the larger 
the number of blades and the smaller the value of y,, rela- 


tive to po. 


Calculation of the Bound Circulation 

The bound circulation is given by the jump in the 
velocity potential at the trailing vortex sheets. Re- 
membering that w/w was taken to be unity, the ex- 


pression for the bound circulation per blade becomes 


z 2m+1 (iH) 


Bol S i, > B 
- — ; a 
27wv ecg. ee ae 
lom4 l (i) Ko», l (f1) 
) B = 3 


Calculation Procedure 

Although the values of a, can be obtained by the 
standard method of successive approximations, the 
labor involved was not considered justified for this 
investigation. Except as a check on the two-bladed 
case an approximate method which will be described 
was used to calculate the circulation distributions pre- 


sented in the figures. 


For the approximate calculation of the change in 
the circulation distribution due to the hub, it is as- 
sumed that, compared to the contribution from the 
b,, terms, the change in the a,, terms can be neglected. 
With this assumption, the bound circulation for the 
finite hub case will be given simply by the sum of the 
b,, terms and the circulation of the zero hub case. If 
A denotes the difference between the circulation with 
and without a hub, then 

Koma+1 . (@) 
Bol 2 > > 8 - 


9 = is Oy > 5 v= 
<TWU 7T m=0 K Im+1 (Mo) 


In order to calculate the b,,’s, consider Eq. (2%). 


b,, is composed of two terms, one involving a,, and the 
0, the mag- 


other involving ii om+1 . (Hn). Form = 


p/») (a) term is about twice to ten 
term while 


nitude of the 7, 
times as large as the magnitude of the a, 


Thus a small 


for larger m, the ratio is even higher. 
error in the value of a,, will not affect the value of 0, 
With this in mind, the approximation 
values 

For 


to any extent. 
is made in the calculation of the b,,’s that the «,, 

remain unchanged with the addition of a hub. 
convenience, these values of a,, given in Goldstein are 


repeated here. 
Qm = —[uo?/(1 + wo?) Am + €n 


where A, = 1, 34, = 1/2, 5Ae = (1-3)/(2-4) 
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(1-3-5) /(2-4-6), and for B ) 
fo = 2; « = —0.061, é, = 0.013, 
co = ¢€3 = os . Q 
My = Be & = 0.047, 1 = 0.007, 
6 =é&=... 0 
bo = 5: é = —0.033, €| 0.004, 
ce = €3 = 0 


and for B > 2 «,, = 0. 

These assumptions considerably simplify the calcu- 
lations. Whether or not they are justified is best 
answered by considering actual calculations which 
have been performed. Circulation distributions have 
been calculated for two-, four-, six-, and eight-bladed 
propellers with hub radii of 0.3 and 0.4 of the propeller 
radius for different values of wp. In addition, for the 
two-bladed propeller, hub radii of 0.1 and 0.2 of the 
propeller radius have also been considered. 

Some idea of the errors involved in the assumptions 
can be obtained from the magnitude of the terms in the 
square brackets of Eq. (30). It is evident that the 
larger the number of blades, the larger uo, the smaller 
u,, or the larger the index m, the closer these terms will 


be to unity. For example, for B = 2, wy = 2, and uy, = 
0.8 the worst variation from unity for any of the four 
bracketed terms for m = O is 5 per cent. Form = | 


this difference is even smaller amounting to only 2 


per cent. For four blades with the same wy and y,, 
the worst difference for m = 0 is reduced to 0.7 per cent 
and for m = 1 amounts to 0.2 per cent. Thus com- 
parable differences in the a,, values would be expected 
between the zero hub case and the finite hub case. 

Another indication of the effect of the approxima- 
tions is the value of the circulation calculated at uw. It 
was found in the case of two blades, that to the accuracy 
of the calculations, the circulation did not go to zero at 
the tip. This, of course, is due to the fact that the 
(»'S Were assumed unchanged by the addition of a hub. 
However, in none of the cases considered did the cal- 
culated circulation at the tip exceed approximately 3 
per cent of the maximum circulation. In such cases a 
first-order correction was made to dp to insure that the 
circulation vanished at the tip. 

The system of Eq. (30) for a,, was solved for dy for 
the case of B = 2, wy = 2, uw, = 0.6 by assuming first 
that aq; = ad = a; = ... = Oand solving one equation 
for da; then assuming ad. = a3; = a, = ...0 and solving 
two equations for a) and so on up to five equations. 


1955 
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The resulting a,'s were expressed as a power series jy 
inverse powers of the number of equations used and th, 
limiting d) obtained as the number of equations becon, 
infinite. In this case a value of a) = —0.863 was ob 
tained as compared to a value of —O.S61 quoted by 
Goldstein for the zero hub case. 

It is rather difficult to say where the approximat 
calculations hold or for what conditions they are in 
sufficient. However, in general the approximat 
method appears satisfactory for more than two blades 
for not too extreme ranges of wy or w,. This means 
that uw, should not be much greater than four-tenths oj 
wo and pw no smaller than one. For the two-bladed 
propeller it appears necessary to solve system (30 
for the a,,’s for uo less than 2.0. 


Results of Calculations 


The effect of a finite hub is clearly evident in Fig 
In this figure circulation distributions are presented for 
a two-bladed propeller operating at uy = 2 having rela 
tive hub radii of one-, two-, three-, and four-tenths of 
the propeller radius. The insertion of a hub is seen 
to have an appreciable influence on the inboard cir- 
culation; particularly at the surface of the hub. For 
example, in the case of Fig. 1, the insertion of a hub 
having a radius of 0.2 the propeller radius increases 
the circulation at the 0.2 radius by a factor of 1.56. 

Comparison of Fig. 2 with Fig. | shows that the 
effect of a hub is not as pronounced for a larger number 
of blades. Fig. 2 presents the calculated circulation 
distributions for an eight-bladed propeller with rela- 
tive hub radii of 0.3 and 0.4 operating at a wy of 2.0. 
In addition to Figs. 1 and 2 results which have been 
calculated to date for two-, four-, six-, and eight-bladed 
propellers are summarized in Figs. 3, 4, 5, and 6. In 
these figures which are for specific values of uo and u,, 
the ratio of the bound circulation with a hub to the 
circulation without a hub is presented as a function of 
the reciprocal of the number of blades for various values 


of relative radii. 
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Simple Shear Flow Past a Fiat Plate in an 
incompressible Fluid of Small Viscosity 


Ting-Yi Li 
Department of Aeronautical Engineering, Rensselaer Polytechnic 
nstitute, Troy, 


March 28, 1955 


\ THE STUDY OF high-speed viscous flow past a two-dimensional 
body it is usually necessary to consider a curved shock wave 
Conse 


tently, there exists an inviscid rotational flow region between 


mitting from the nose or leading edge of the body 





the shock wave and the boundary layer. Such a situation arises, 
for instance, in the study of the hypersonic viscous flow past a 
flat plate. The situation is somewhat different from Prandtl’s 


classical boundary-layer problem. In Prandtl’s original problem 
the inviscid free stream outside the boundary layer is irrotational 
while in a hypersonic boundary-layer problem the inviscid free 
The possible effects of 
vorticity have Libby 

In the present paper, the simple shear flow past a flat plate in a 


stream must be considered as rotational 


been recently discussed by Ferri and 


fluid of small viscosity is investigated. It can be shown that this 


problem can again be treated by the boundary-layer approxi 
nation, the only novel feature being that the free stream has a 


constant vorticity. The discussion here is restricted to two 





limensional incompressible steady flow. The basic equations 
ind boundary conditions are expressed in the usual notations as 
follows 
u(Ou/Ox) + v(0u/Ov) = v(02u/dy? (1 
(O0u/Odx) + (dv/dy) = 0 (2 
y = QO, u=v=0 (3 
yo &, u = Uo ou/oy = +w (4) 
where 
u = uo X= woy (5 
Eq. (5) defines a simple shear flow in the undisturbed free stream, 


being constants. The vorticity in the free stream is 


is satisfied if a stream function y¥(x, y) is 


M ind Ww 
equal to w Eq. (2 
introduced as follows: 


u = Oy/oy —v = Oy/Oox (6 


Then the boundary-layer equation and its boundary conditions 


become 
OY Ov) (07y/Oxdy) — (Ov /dx) (07y/dy?) = v(O8Y/dy3) (7 
y=(0 Oy/Ox = OY/oyvy = 0 (Ss 
= Oy /Ov = u 0*y /Oy? = +w (9) 
Now rewrite the free-stream velocity distribution as follows: 
ul = My [Ll (wy/ty/V vx/uo) (y/V vx /uo)] (5a) 


lhere appear two nondimensional parameters in Eq. (5a 


n = (YV o)/V vx Blasius variable (10) 
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oO} papers pub 
Publication is completed 


hold itself 


doe § not 


E = wo/Uy/V vx/u ‘vorticity’ number 11 


¢ 


While 7 is a 
a new parameter in the present problem 


familiar variable in boundary-layer studies, & is 


Since as y > ©, the 
present boundary-layer flow must join smoothly the free stream 
as given in Eq. (5a), it seems natural to assume that the solution 


must depend on the ‘‘vorticity’’ number in some manner. In 


deed, uc/V vx/uo ~ u/d5, where 6 denotes a boundary-layer 


thickness, represents essentially the average vorticity strength 


inside the boundary layer. Therefore — can be interpreted as a 


ratio of the free-stream vorticity and the average vorticity in the 


boundary layer. The boundary-layer approximation implies 


that a thin layer of strong vorticity exists adjacent to the plate 
surface. Hence the ‘‘vorticity’’ number must be small—i. e 
gé< 1 The stream function ¥(£, 7) can therefore be expanded 


as the following series: 


Y = vV wWox/v lfiln) + Efi(m) 4 (12 


» 


Substitution of Eq. (12) into Eq. (7) yields 


uo?/x[—(1/2)fofo” + (1/201 2 
(u xX) Uf rel + (13 


where d/dn is denoted by the prime and terms of 0(£?) and higher 


order terms are indicated by the dots. Equating like power 
terms in from Eq. (13) vields 
ff." + OF! =O (14 
2f,'"! + fof” — fo'fr’ + Bhe’h = 0 (15 
The boundary conditions become 
7 =U i f = () = f = () (16 
oe fy = 1 f° = 1 17 


It is noted that fy is the Blasius function and that Eq. (15) is 


linear in f;. The solution for f; can be obtained by the usual 
technique. First, taking f, © (1/2)an? for 0 < » < m, one finds 


the solution f; satisfying the conditions of Eq. (16 


, ; Van3/3 ” 
fy = (2/3)Af f ‘ T2/(¢)de/¢ a = 0.332 


. 


as follows 


first kind of 


~ 


is the modified Bessel function of the 
©, taking f 


where /2/,(¢ 
the order 2/3. On the other hand, for 7; < » < 
n — 8, one finds the solution f; satisfying the conditions of Eq 
(17) as follows: 


P a ] n — B)* 
"age | 2B | , n — B)?/4] F 1, -, ) dn — 
n 2 4 


i, = 


where , F(a, c, x) denotes the confluent hypergeometric function 


The integration constants A, B, and C can be determined uniquely 


by requiring that fi, fi’, and ff,” evaluated at » = m from both 
Eq. (18) and Eq. (19) are equal. The skin friction 7+ can be 
expressed as follows: 

rT = p(Ou/Oy)y=0 = (pmo? ?/V wx fy"(O)  Ef,"(0)] (20 


It follows that the local skin-friction coefficient can be given as 
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Fic. 1. The undisturbed velocity in the free stream (/eft, type 
I; right, type II) 


g 
Cy = = 2V v/uox [fo"(O) + Efi"(0)] (21) 
(1/2) puto? 
In this last formula, fo’(0) = @ = 0.332, f;"(0) can be found from 


Eq. (18) as 


P sats 
fh"(O) = - | (22) 
$(3 )/381(2/3) 
Taking 7, = 3.0 and using the results of Eqs. (18) and (19), one 
finds that A = 3.28. Thus Eq. (22) yields f,"(0) = 0.278 
Hence, Cy becomes 
Cy, = (0.332 + 0.278¢)2V v/uox (21a) 


Therefore a shear flow of the types I or II (Fig. 1) will tend 
respectively to cause an increase or decrease of the skin friction 
on a flat plate. Suppose that in a certain region on the plate 
surface the average vorticity inside the boundary layer is, say, 
ten times stronger than the free-stream vorticity—i.e., — = 
0.1—then the local C; will be affected by the free stream vorticity 
about +8.4 per cent. The displacement thickness can also be 
calculated as usual. The vorticity effects are again a contribu- 
tion of a term of the order of &. If the present computations 
are to be carried to £? terms, it is necessary to consider the pres- 
sure change induced by the displacement effect of the boundary 
layer. Such a refinement, however, has not been carried out at 
the present. It remains to say a few words of the limitations of 
the present approach, viz., the series expansion method will be 
However, it is doubtful whether the 


inaccurate for & = 1. 
boundary-layer approximation is still applicable when the aver- 
age vorticity in the boundary layer and the free-stream vorticity 


are of comparable order of magnitude. 
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Relative Size of Turbojet Engine Components 


D. Willens* 

Aerodynamic Engineer, Propulsion Research Corporation, Santa 
Monica, Calif. 

March 31, 1955 


O™ OF THE IMPORTANT FACTORS to be considered in turbojet 
engine design is the engine frontal area. As flight speeds 
are pushed farther into the supersonic regime, the effect of engine 


* Now, Head, Aero-Thermodynamics Group, Turbo Research Corpo 


ration, Santa Monica, Calif 


SEPTEMBER, 1955 


frontal area increases in importance. Thus, it becomes neces 
sary to determine which components are most likely to fix th 
engine outer diameter. In addition, it is necessary to deter 
mine the effect of engine design flight Mach Number upon th, 
component relative sizes 

A brief analysis of this problem was presented in reference 
for nonafterburning engines. In this note the analysis has bee 
based upon a different set of initial assumptions and extended t 
include the afterburner and discharge area of the exhaust nozz\ 

From the continuity equation, the flow area at any station j 
the engine can be determined in terms of the area at some refer. 
ence station from the ratio of the densities and flow velocities 
the two stations. The reference area used in this paper is th 
compressor inlet frontal area (see Fig. | 

The density at each station in the engine can be determine 


from the thermodynamic cycle as a function of altitude, flight 


Mach Number, compressor pressure ratio, turbine inlet tem- 


perature and the exhaust nozzle temperature after reheat. Th 


flow velocities, and in some cases the radius ratio, must be stipu- 


lated in order to obtain the frontal area of each component 


Although the results presented here are for the engine design 


point, consideration must be given to the off design operation i: 
order to draw generalized curves which will be applicable over a 


wide Mach Number range. It has been assumed that the en- 


gines resulting from the design point analysis would be operated 
with constant mechanical r.p.m. at off design flight conditions 
Therefore, it was necessary to allow the assumed value for the 
per cent of sonic flow at Station 1 to vary with flight Mach 
Number in order to prevent a compressor inlet choked condi 
tion at an off design operating point 
variation in compressor inlet per cent of sonic flow with flight 
Mach Number is shown in Fig. 2 and is based on a value of 6 
per cent at sea level static. The remainder of the assumptions 
used to determine the frontal area of the combustor, turbine 
afterburner, and exhaust nozzle discharge are listed in Table | 
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Fic. 1. Schematic engine diagram showing location of com 
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TABLE | 


A) Compressor 


1) per cent sonic flow in frontal area Fig. 2 
B) Combustion Section 
1) full annular combustor 
average inlet velocity, ft. per sec 100 
8) jnner-to-outer radius ratio 0.4 


C) Turbine 
1) average total temperature at first stage nozzle 2,060°R 
& 


2) first stage nozzle hub ratio 


3) first stage nozzle discharge angle 25 
1) first stage nozzle choked 
5) last stage discharge axial Mach Number 0.5 
6) constant hub diameter 
D) Afterburner 
1) average inlet velocity, ft. per sec $50 
E) Exhaust Nozzle 
1) average total temperature 3,500°R 


complete expansion to ambient pressure 


In actual practice these quantities may be varied somewhat 
in order to obtain an optimum design for a given flight applica- 
tion. However, in order to show the important trends due to 
flight Mach Number and compressor pressure ratio, the values 
in Table 1 have been held constant 
5, and 6 show the effect of compressor pressure ratio 


Figs. 3, 4, ° 
PR.) and design flight Mach Number on the ratio of frontal area 
of the various engine components to the frontal area of the com- 
pressor. It should be noted that the compressor pressure ratios 
shown are the values which exist at the various flight Mach 
Numbers 

Consider first the subsonic flight condition of 1J = 0.9. For 


i PR, = 6, the turbine is the largest component. Increasing 
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Fic. 3. The ratio of combustor frontal area (Station 3) to 
compressor frontal area 
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Fic. 4. The ratio of turbine frontal area (Station 5) to com- 
pressor frontal area 
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compressor pressure ratio decreases the size of all components 
At PR, = 


12 the compressor is slightly larger than the remaining com 


relative to the compressor except the exhaust nozzle 


ponents. 

Consider secondly a flight Mach Number of 2.0. At PR, = 
4 the largest component is the completely expanded exhaust 
nozzle. Increasing compressor pressure ratio to & decreases 
all component sizes except the exhaust nozzle discharge 


area 

Consider finally a flight Mach Number of 3.0. At PR, = 2 
the combustor is the largest component, being 2.24 times the 
compressor frontal area. This component can be reduced in 
size by designing the combustor for an inlet velocity greater 
than the initially assumed value of 100 ft. per sec. Using a 
higher combustor velocity, it is seen that the exhaust nozzle is 
now the largest component at 1.92 times the compressor area 
At PR. = 4, the largest component is the exhaust nozzle dis 
charge 

To summarize the above results, the frontal area of high pres- 
sure ratio turbojet engines designed for an altitude flight Mach 
Number of 0.9 will be determined by the compressor. For 
afterburning engines designed for flight Mach Numbers from 
2.0 to 3.0 and using an exhaust nozzle with complete expansion 
to ambient pressure, the maximum frontal area will be at the 


nozzle discharge. 
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Application of Kelly’s Theory to Cone-Cylinder- 
Frustum Bodies of Revolution* 


Robert E. Lavendert 
Redstone Arsenal, Huntsville, Ala 
March 9, 1955 


A‘ IMPROVEMENT IN THE THEORY of viscous crossflow over 
bodies of revolution has recently been made by Kelly.' 
Excellent correlation of the normal force and pitching moment 
variation with angle of attack was obtained between experiment 
and theory for a series of bodies in both subsonic and supersonic 
flow. The theory assumes that the axial distribution of the cross 
flow drag for an inclined body is similar to the drag for a circular 
cylinder moving crosswise in a fluid when started impulsively from 
rest. In addition, the proper cross-flow drag coefficient is as 
sumed to depend upon whether or not the axial flow is laminar or 
turbulent. 

Kelly has applied the theory to cone-cylinders and ogive 
cylinders for cases of essentially all laminar or all turbulent axial 


“ Statements and opinions advanced in this note are those of the author 
and not those of Redstone Arsenal 

7 Acting Chief, Aerodynamics Branch, Technical Feasibility Studies 
Office, Ordnance Missile Laboratories 


SACy 


1 
: = sin? a [ rCp dx 
o J0 


and 
SdAC y Jee Pe : 
- : = sin? a rCpx dx 
7 J 0 

where 

ACy = additional normal force coefficient 

ACy = additional pitching moment coefficient about nose 

S = reference area 

d = reference length 


The cross-flow drag coefficient Cp as presented by Kelly is 
. . x 
Cp = Chg A tana +B 

; r 


The steady state value Cy, is considered to be 1.20 for axial laminar 


(1) ForO <x <Sa;r =x tand,n 0.5 


Pree ; I =F Aa? $Ba* 
sin? @ rtp dx = Chg = Aa?a’ + 2Baita'® 4 T es 


(2) Fora<x<br=n 0.5 


r 
‘ ‘ B : 
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b l 
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Ja i 
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Cone-cylinder-frustum body of revolution 


flow. The application of Kelly’s theory to cone-cylinder 
frustum bodies (see Fig. 1) is presented in this note and may bs 
used when the axial flow is all laminar, all turbulent, in transition 
near the cone-cylinder juncture, or in transition near the cylinder- 
frustum juncture. 

The expressions for the additional normal force and pitching 
moment due to crossflow may be written 


b ~, 
: rCp dx +4 } rCp dx 
Ja Jb é 
2h » 
+ } rCpx dx + } rCpx dx 2 
ae Jb 


tan?’ a + D tan® a 3 


yd 


flow or 0.35 for axial turbulent flow. 
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Eqs. (1A), (1B), and (1C) are evaluated separately so that the properly Cp, can be used in each region Phen, for any given body, the 
three equations may be combined in Eq. (1) to obtain a single equation for the addtional normal force In obtaining the expressions as 
. function of angle of attack (radians), use was made of the expansions 
g sin? a tana = a + (a?/15) 4 
—j sii- a tan’ a = a + 2a'/3 t 


sin? a tana = a’ 





| The evaluation of Eq. (2) is made in a similar manner 
$) ForO <x Sa;r =xtand,r 0.5 
wae { Aa SBa 16Da 
— sin? a rCpx dx = Cp, Aaa +4 Baad 4 : : a 2A 
' - 3 3 15 q 3 


For « x<b;r=r O05 
‘i i : jl } A(b* — as SB(b> — a 16D | 
I~ @ rCpx d = ( De : 1(b° — a®ja _s Bib’ — a®)a ¢ 7 + = + - a 2B 
linder : 13 5 L 45 15 J 4 
lay be 6) For o:% , x — 5) tan w,? O.5 
sition 
linder- ‘ ' jl ic? — 0 2 7) 
sin? @ rCyx dx = Cp, \-= A(c? — b®)a* + Blza® 4 = + BI; + DI; | a 2C 
It 13 15 3 : 
tching 





The coefticients of a’, a®, and a’ in Eqs. (2A), (2B), and (2C c’/c = 1 Cri 3) 
re also evaluated using the proper Cp, for each region. The A] ’ ee 
: 3 ; : ; Also, ¢ = ¢c,andc, = ( 
quations may then be substituted into Eq. (2) so that a single Using conventional NACA shorthand notation,? the lifting 
(3 equation is obtained for the additional moment coefficient “ae :; ; ; : 
effectiveness of a leading-edge flap may be expressed in terms 


Use of the above equations enable one to correlate experimental ; . es ; : 
of linearized trailing-edge flap parameters as follows 


lata with the theory at low Reynolds Numbers and then to extend 
the calculations to higher Reynolds Numbers ci = Cia ago 
= ( a’ + §'(1 Ty a 
REFERENCE 
Differentiating with respect to 6’ at constant section lift 








(1A ' Kelly, Howard, R., The Estimation of Normal Force and Pitching Mo 
ent Coeficients for Blunt-Based Bodies of Revolution at Large Angles of (0a’/06")-. = a = l a 4) 
Attack, NOTS Tech. Memo 998, May, 1953 ; 
Differentiating with respect to 6’ at constant a 
(1B = + Cl, — as Cl, 4 
In a like manner, it can be shown that at constant section lift: 
Estimation of the Lift and Moment at ag sa 6 
Parameters of Leading-Edge Flaps 
(1c 
Malcolm J. Abzug 
Design Specialist, Aerodynamics, Douglas Aircraft Company, AIRFOIL WITH MTOR. WIN 
Inc., FI Segundo, Calif TRAILING-EDGE FLAP LEADING-EDGE FLAP 
April 11, 1955 
© 
7 Pee AERODYNAMIC LIFT and moment parameters of leading- c™— § ‘(SHOWN NEGATIVE 
+) edge flaps may be obtained from trailing-edge flap theory! or bee 
Irom tests of wings with trailing-edge flaps by means of a simple 
geometrical transformation. This transformation is illustrated 
in Fig. 1, where two flapped airfoils presenting identical profiles 
to the relative wind are shown. Letting primed values refer 
to the leading-edge flap case: WIND > 
Oo 
a=arté (1) . : ; ; , ; . . , 
Fic. 1. Flapped airfoils presenting identical profiles to relative 


6° = —6 (2) wind 
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Using Eqs. (3) and (4), the leading-edge flap lifting effectiveness 
parameter as’ is compared with the theoretical trailing-edge flap 
value for incompressible flow’ in Fig. 2. This figure shows that 
the lifting or control effectiveness of narrow-chord leading-edge 
flaps at subsonic speeds is many times less than that of narrow- 
chord trailing-edge flaps 

Comparisons of theoretical and experimental section lift and 
moment coefficients for two leading-edge flaps’ ® are presented 
in Fig. 3. Reasonably good agreement is indicated. 
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Remarks on ‘‘An Exact Solution for Two- 
Dimensional Linear Panel Flutter at 
Supersonic Speeds’”’ 


S. F. Shen 
Associate Professor, Department of Aeronautical Engineering 
University of Maryland, College Park, Md 


April 25, 1955 


r A RECENT NOTE,! Goland and Luke gave an ‘‘exact solution’ 
of the panel flutter problem at supersonic speeds, based upon 
the integro-differential equation of motion employing the linear 
ized aerodynamic pressure without other assumptions. The 
solution in Laplace transform was neatly formulated, although 
However, certain conclusions 
Because of the ‘‘exact”’ 


numerical results were not given 
were drawn near the end of reference 1. 
nature of the treatment, these conclusions are likely to receive 
For this reason, it is felt that some clarifying 
We shall show that their conclusions ar‘ 
A more 


wide attention. 
remarks are in order. 
open to question by means of a simplified example 
careful attitude toward assessing flutter instability will also be 
emphasized. 

In reference 1, Goland and Luke defined a parameter 7 which 
is roughly proportional to the ratio of the dynamic pressure to the 
inertia running load at the fundamental natural frequency. For 
y>>1, the important conclusion was that flutter of both plate 

“Tt can also be shown without difficulty 


and 


membrane cannot occur. 
that flutter cannot occur for the membrane if the Mach Number 


M is sufficiently large.’”’* It should be noted here that they were 


* Quotations in this note are all taken from reference | 
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secking the mode shape in which the panel is capable of a 
By ‘‘flutter’’ was meant 


reall 
self-excited simple harmonic motion 
the existence of such a mode 

In this light, then, the nonexistence of a mode shape capable of 
<elf-excited simple harmonic motion, or at least an almost simple 
harmonic motion, at any y, is indeed physically difficult to visual 
2 By going to higher and higher modes, hence frequencies, the 
inertia loading seems always able to attain such a magnitude as to 


x of the same order as the elastic and aerodynamic loadings 
And, for the very highest modes, the inertia loading can always so 
ompletely dominate the aerodynamic loading that the normal 
modes and frequencies at the high end of the spectrum would not 
ye seriously affected by the presence of air More precisely, we 
mean here that the asymptotic behavior for the higher modes at 


riven y Should be [w(y)} /@(O)] ~ 1 + o( 1), where w(y) is the 
frequency for the given y. The term 0(1) may of course be com 
plex, indicating some amplification or decay with time. But a 


mode would still exist The conclusion of reference | is appar 
ently contradictory to this reasoning 
the explanation probably lies in the limiting 


Analytically, 
While not stated in 


process which led to Eq. (19) of reference 1 
letail, the derivation of Eq. (19) must have been made by for 
in Eq. (14) but keeping all other parameters 


mally letting y > 
as the eigenvalue) constant. Such a 


including the frequency 
rules out the possibility that the frequency 


proc edure ck il ly 
however 


parameter may grow without limit, given a fixed y, 


large it may be 
lo illustrate the inadequacy of such a limiting process, we 


onsider the following elementary problem: 


—_ y2’ + g2p2- = ()/ , 
{ ) 
(0) = 0, 2(2) = ON 
vhere symbols in reference 1 are used. This corresponds to 


essentially a certain membrane problem with quasi-steady aero 


namic forces. Applying the Laplace transform, 


(p? — yp + a*k?)g = s_, (2) 
or 
5(p)/t,. = (Pp? — vp + ak?) (3) 
Letting formally y ~ o, 
(p)/2, = Pp — y)™ (4) 
or, finding the inverse transform, 
(x)/2,, = —(1 — e7*)/4 (5) 


Comparing against Eq. (19) of reference 1, we see that again here 
the dominant term is e?7*/y. One would have to conclude like 
wise that ‘for large y a physically feasible flutter mode is non 
existent.”’ 

Actually, this elementary problem can be solved explicitly to 


vield 


v)/s = (2/nmr)ev"/- sin (nwx/2) om 1, 2.3, , (6 


with 2ak = Vo n?x? + y?, to satisfy both end conditions. The 


fact that k increases with y reveals itself. It may be noted that 
our elementary problem belongs to the classical Sturm-Liouville 
problem, for which the behaviors of cigenvalues and functions 
ire well-known. The direct calculation above is merely to bring 
out the result in a transparent manner 

At the risk of causing a controversy, we take this opportunity to 
emphasize that, even if a possible mode should exist for the 
issumed simple-harmonic motion, without other clues it may not 
be correct to conclude that there is flutter instability. Obviousl) 
here only the condition for marginal stability is fulfilled. Our 
elementary problem defined by Eq. (1) and solved by Eq. (6) is 
an example In this example, regardless of the value of y, its 
flect, provided re il, is simply to modify the mode shape and the 
frequency which goes with the mode. There is a continuous 
change of the mode shape and frequency with a continuous 
change of the parameter y from zero upward. The motion is 
ilways neutrally stable. It corresponds to nothing more than 


{ie natural vibrations with y + 0 
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fo contirm flutter instabilitv, our belief 1s that divergent 
motion must be shown to exist on one side of the “flutter bound 
ary.”’ The job is admittedly not trivial, but probably essential 
in panel flutter since our example suggests that wrong interpreta 
tion may result otherwise We have made certain investigations 
along this line, proposing to treat panel flutter as a complex 
eigenvalue problem.* It is hoped that the results will be sum 
marized in a forthcoming report 
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On the Flow Behind an Attached Curved Shock 


S. |. Pai 
Associate Research Professor, Institute for Fluid Dynamics and 
Applied Mathematics, University of Maryland, College Park, Md 


April 28, 1955 


IT A RECENT PAPER,! Eggers, Savin, and Syverston mention 
that in my analysis? of the flow behind an attached curved 
shock, the entropy gradient is neglected. This statement is defi 
nitely untrue. In my final results, [Eqs. (59) and (60) of refer 
ence 2] the effect of entropy variation is included 

The essential feature of this problem is shown in Fig. 1. The 
disturbance at any point, say ‘‘3"’, may be roughly decomposed 
One is the basic disturbance which depends 
and the other 


into two parts: 
mainly on the local slope at 3—i.e., 3-4 in Fig. 1 
is the reflected disturbance—i.e., 2.3 in Fig. 1—which depends 
on the so-called ‘“‘reflection factor’’ at the shock (point 2) as well 
as the basic disturbance from 1—i.e., 1-2 in Fig. 1 Eggins, 
Savin, and Syvertson try to show that the reflection factor is 
small so that the reflected disturbance 2.3 is negligible in com- 
parison to basic disturbance 3-4. But they forget that the re 
flected disturbance depends not only on the reflection factor but 
also on the disturbance 1-2. If the disturbance 1-2 is much 
larger than 3-4, the reflected disturbance 2-3 may not be small 
in comparison with 3 4 even though the reflection factor is small 

In my analysis, I did not separate the disturbance locally in 
the manner shown in Fig. 1, but considered simultaneously the 
whole shape of the body and the entropy variation. In Eq. (59), 
the reflected disturbance includes both the effect of the so-called 


“reflection factor’’ of Eggins, etc., and also the incident disturb 





SHOCK 





FLOW BEHIND AN ATTACHED 
CURVED SHOCK. 


Fig. I: 



































658 JOURNAL OF THE AERONAUTICAL SCIENCES SEPTEMBER, 1955 
LNCes In this manner, we can not and should not use a singh iis 
reflection factor to determine the magnitude of the reflected dis | : t 
turbance M=| p ; p 
I had discussed this point with Mr. Syvertson in several com n —t (1) . 4 5 > 2 ; 
munications. My “reflection factor’? is not the same thing as 1 M5 ” 
his “reflection factor." My reflection factor is so defined that ‘ fi eee a Le _ ‘ TTED 
it is convenient to use in my analysis. I do not believe that any a WAL 
single reflection factor could be used to determine the total re ite ; 
flected disturbance One should also consider the effect of in a 6 " o ai 
cident disturbance simultaneously —i.e., the shape of the body in FIG | 
front of the point considered. Eggins, Savin, and Syvertson’s 
inalysis can at most apply to bodies whose slopes at every point — — 4 
are almost the same—i.e., disturbance 1-2 is almost the same as Vs " v, 
3-4. In my analysis, there is no such restriction and I did not . y ’ Y ‘ Y . 
“partially neglect’? the entropy gradient In some numerical oe w | 
examples,* I found that the resultant disturbance differs con 
siderably from the basic disturbance which shows that in those . i bs t+ ; . I 
cases, the reflected disturbances are not negligible 
REFERENCES 
tee 
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+ 
j - 
Free Streamline Analysis for Slotted 
Two-Dimensional Supersonic Minimum- ~- Stil 
Length Nozzles 
Robert Wesley Truitt 
Professor of Aeronautical Engineering, Virginia Polytechnic 
Institute, Blacksburg, Virginia Fro 
May 10, 1955 
nd ire 
ce ATTAINMENT of a given supersonic Mach Number in a \ 
constant cross-section two-dimensional slotted test section p 
would require, under minimum-length conditions, a certain ratio : V 7 Using 
of slot width to working section width. The condition for the p 


absolute minimum-length symmetric nozzle is shown schemati 


cally in Fig. 1. The angle v is the Prandtl-Mever deflection such " 
that kK | » Substit 


y p’/2 (1 


vy’ being the deflection of a sonic stream necessary to produce a F | G ? b 
. 


given working section Mach Number J/’. The solid lines show 
the conditions of the solid-wall supersonic nozzle; the broken 





line represents the transonic-type constant cross section with 
fixed-width slot. The length, d, is the minimum possible for 17’ 6 es Pe Fes - 





It is clear from Fig. 1 that, if the nozzle section is kept constant, 
a Mach Number, J/’ > 1, would require a certain mass rate of 





flow through the slotted wall. The fluid mass will have the , 16 
velocity, v, normal to the slotted wall 4 Sa Tu 
2), 18 
SINGLE OR MULTISLO1 L/w 1 Pr 
: : ; id ; : s show 
A cross section of the ideal flow conditions for a multislotted : 
; : f eities 4 Imiting 
wall in region 2 may be depicted as shown in Fig. 2a Due to _ 
nanan siphon 2 CRITICAL to 
symmetry, the flow may be simplified to that in Fig. 2b, where z , , 
‘ ° ° e » 1 Teache 
is taken as the tunnel width; Z is the actual slot width; L’/L J Ms | 87 
n, the contraction coefficient / 
The analysis will be made assuming the flow through the slot (s) A 
O na 4 1 4 —a B 





is the same type as the compressible free streamline flow through a 


| 10 2 i4 mi 16 i8 20 


In a transonic-type working section, the pressure, p’, and 
density, p’, in the pressure chamber are the same as in the test F | G. 3 


two-dimensional orifice treated by Busemann 
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ne; clio fore, for any cross section of the slotted wall in applied to get 
t ¢ may apply the continuity equation as (see Fig 
p V ve? + [Qy/y — 1] }(pe/s } 
pool pu L 
Eq $) assumes isentropic flow between region 2 and _ the 
TED OF = pwvatl p'v'nl - pressure chamber so that 
ALL : 
her ire the mass density and normal velocity com 
1 tively, in region 2 f jl y — 1)/2].M"\ 
rl | ratio of tunnel width to slot width is then D Bl vy — 1)/2)Mo24 
w/L = nl p'/p2)(v' /2 3 oa 
— A enterline stream line, between region 2 and the 
ber, Bernoulli's compressible equation may be po /f »» / p (j 


Using kq 1), (5), and (6), Eq. (8) becomes 


u fl + [(y — 1)/2] 12}! 2yp fl + [(4 2). , 
= 9 - 1 + : | ‘ 
Lt (1 + [(y — 1)/2}."24 \ y — 1)es2p’ (1 + [(y — 1)/2) { 


iction coeflicient for the two-dimensional orifice may be written for the present case as 


olution for the compressible contr: 





From Eq 1), (5), and (6) we may write 
c’ ( y—1 1+ [(y -— 1)/2|M { 


Using Eqs (8) and (9) and the relation for the speed of sound at L’ as 


i 


! Kg pecomes 





v2 \? 2 [+ iy = : 
r+2 1- (*) - — | 
\ ( y—1 Q1+ [(y — 1)/2] M2? f 


V2 sin »y = MnV yp2/p2 sin v 12 
id from Eq. (10 
c’ = VV aY (p p' (p’ f Sin pv 14 
sing Eqs. (5) and (6), Eq. (18) may be written as 
; be Mii aye 
Vo | Y 1)/2)].M 14 


V2(sin v 
c "VSG - DAM: 


14) into Eq. (11) gives the required relation 


Substituting Eq 
pl + [Cy — 1)/2).Me7{! 2 lj 1+ [(y — 1)/2] Me?) 
a i 4 > a. - 
(1 + y — 1)/2|."24 y¥ - I \sin* » Wo? | Ll + y — 1)/2|M'24 - 
! 2 ¥ $l t ly - DAI") 2 §1t+ ((y — 1)/2|M j 
1 + 1 — (sin? vy) 1? _ = - | 
ae 1 1 + [(y — 1)/2) 174 y—1 (14+ [ly -— 1)/2)M f 
—t 
Eq. (15) is the ratio of the working section width to the slot width . 
isa function of AJ’, 2, and vy. The Mach Number, JJ, in Eq : 
15), is the Mach Number the sonic stream will attain in region 2 Buckling Efficiencies of Plate Materials at 
t Prandtl-Meyer expansion of vy = v’/2. A plot of Eq. (15 Elevated Temperatures 
sshown in Fig. 3. It is interesting to note that Eq. (15) has a 
limiting val va - 2 9 Th; ws , ). W. Clark 
Imiting value at J 1.87 (see Fig. 3 Phis critical value is ‘sega! % P 
- 1 . ’ ° . 5 rcn cCngir er tCngineering Vesigr U r UIT 
‘ue to the fact that the Mach Number (v’/c’) in Eq. (8) has Researc . - ag : Pag 
reached th : : Research Laboratories uminum Company of America 
ichy 1e sonic value. New Kensington, Pa 
April 30,1955 
REFERENCI 
20 Busemann, A., Hodographenmethode der Gasdynamik, ZAMM. Vol. 17 ea. rHE JOURNAL, Gerard! compared the structura 
1 efficiencies of plate materials at various temperatures. He 
leveloped two structural efficiency parameters,” one that 


ipplied to plates that fail because of elasti buckling and the 
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other to plates that buckle at the compressive yield stress. The 
purpose of this note is to show that the latter parameter does 


not give a true comparison of structural efficiency 


In comparing the weights of two plate structures of different 
materials which fail by plastic buckling, Gerard assumes that 
the load and the dimensions of the plate panels are the same for 
the two cases. The thicknesses of the two plates are adjusted 
so that each buckles at the compressive yield strength of the 
material. The weight is then proportional to the parameter 
p/E,¢*, in which p is density and E£, is the secant modulus at 
the compressive yield strength 

The error in the above reasoning lies in the assumption that 
one variable, the thickness, can be adjusted to achieve two 
simultaneous results: (1) both plates buckle under the same 
load, and (2) both plates buckle at the compressive yield strength 
In general, if both plates buckle at the same load, the buckling 
stress may coincide with the compressive yield stress for one 


material, but it will not for the other 


If two plates were designed so each would fail at the compres 
sive yield strength, the relative weights of the two plates would 
be directly proportional to the density p and inversely propor- 
tional to the compressive yield strength. Thus the weight 
strength efficiency for plates which fail at the compressive yield 


stress can be measured by the parameter, p/a,, 


However, the parameter p/o,, cannot serve as a measure of 
compressive buckling efficiency for the same reason that the 
parameter used by Gerard is not valid. That is, one cannot 
design plates of two different materials to buckle at the com 
pressive yield strength under the same load without varying the 
panel dimensions as well as the thickness 

It may be that no single parameter gives an adequate picture 
of the relative efficiency of various materials in resisting buckling 
in the plastic stress range, since the comparison will depend to a 
large extent on the conditions encountered in an individual de 


sign problem 
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A Note on ‘‘A Study of Surge in Compressors 
and Jet Engines’’' 


W. C. Randels 
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T" reference 1, a mathematical analysis of compressor surge is 
This analysis is probably as complete as is prac 


presented 
tical, but the resulting criteria are so complicated that it is diffi 
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Errata—Nonlinear Lift and Vortex Separation 
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June 30, 1955 


REFERENCE 1 


ere 213, Eq. (3): Change the symbol ‘‘e’’ on the left-hand 
side to ‘‘e”’ (the conjugate of e). First line after Eq.(5): The 
phrase ‘‘By virtue of the symmetry in the flow pattern’ is in. 
appropriate; omit. Eq. (6): Change the variable (2at/a)?4 
to (2aUt/a)?/3 

Page 214, second line in the first paragraph of Section (3); 
The word “‘infinity’’ should be replaced by the phrase ‘‘in the 
streamwise component.’” Eq. (9): Change the symbol “P” og 
the left-hand side to “I” (the circulation). Eq. (10): Change 
all symbols ‘‘a’’ on the right-hand side to ‘‘ax.”’ 


REFERENCE 2 


Page 222, Eq. (4.4): The equation should be written as 
_aP > dt 


"tb (0 = »P Pg 
; dx ‘ dx 


Page 2238, Eq. (5.1): A superscript minus ‘‘—"’ should be added 
to the symbol ‘‘y”’ to denote the value of the perturbation potential 
on the compression side of the wing. Eq. (5.2): <A bracket sign 
‘*)”’ should be added at the end of the equation 

Page 224, Fig. 5: The exponential power 5/3 appeared in the 
formula of C,* should be written as 2/5 

Page 225, the sixth line from the beginning of Section (7): 
The numbers designating the equations are not correct. The 
correction reads, ‘Eqs. (4.8), (4.9), (6.2), and (6.3).”’ 

All symbols ‘‘e’’ appearing after Eq. (2.11), page 220, are mis- 
prints of ‘‘s"’; refer Figure 2, page 219, for clarification 
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cult to see the physical meaning. A more drastic simplification 
which still applies to most actual cases of surge is presented in 
reference 2. This simplification also applies to other oscillator 
problems. The reasoning ef reference 2 was presented earlier in 


reference 3 but this is a little difficult to get hold of 
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